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We study the elastic behaviour of a supercoiled DNA molecule. The simplest model is that 
of a rod like chain, involving two elastic constants, the bending and the twist rigidities. Writing 
this model in terms of Euler angles, we show that the corresponding Hamiltonian is singular and 
needs a small distance cutoff, which is a natural length scale giving the limit of validity of the 
model, of the order of the double helix pitch. The rod like chain in presence of the cutoff is able to 
reproduce quantitatively the experimentally observed effects of supercoiling on the elongation-force 
characteristics, in the small supercoiling regime. An exact solution of the model, using both transfer 
matrix techniques and its mapping to a quantum mechanics problem, allows to extract, from the 
experimental data, the value of the twist rigidity. We also analyse the variation of the torque and 
the writhe to twist ratio versus supercoiling, showing analytically the existence of a rather sharp 
crossover regime which can be related to the excitation of plectoneme-like structures. Finally we 
study the extension fluctuations of a stretched and supercoiled DNA molecule, both at fixed torque 
and at fixed supercoiling angle, and we compare the theoretical predictions to some preliminary 
experimental data. 

PACS numbers: 87.15. By, 61.41.+e 



I. INTRODUCTION 



During the last few years. Single Molecule Biophysics has become a very active field of research. Among the new 
explored areas, one finds the study of elastic properties of biopolymers, in particular the DNA molecule, under physical 
conditions close to that encountered in living organisms |^-|^. The theoretical work reported in the present paper 
was motivated by the results of micromanipulation experiments being performed at Ecole Normale Superieure |^ 
with an isolated DNA molecule, immersed in a solution of given salinity at room temperature. One extremity of the 
molecule is biochemically bound at multiple sites to a treated glass cover slip. The other end is similarly attached 
to a paramagnetic bead with a radius of a few microns which is performing a Brownian motion in the solution. An 
appropriate magnetic device is able to pull and rotate at the same time the magnetic bead. Because of the multiple 
attachments of the molecule extremities, the rotation of the bead results in a supercoiling constraint for the DNA 
isolated molecule. 

When no rotation is applied on the bead, the force versus extension curves are very well described Q by a simple 
elastic model, the so-called Worm Like Chain (WLC) model which describes a phantom chain using a single elastic 
constant, the bending rigidity This model reproduces the data within experimental accuracy, in a very wide range 
of pulling force from F ~ 0.01 pA^ to _F ~ lOpN. One has to go up to 70pN to leave the elastic regime and see a 
sharp increase of the molecule length which is associated with a transition to a new molecular phase, the "S-DNA" 

@- . . . 

In contrast, in stretching experiments performed with supercoiled DNA molecules 1^,^, two molecular structural 

transitions appear within a more restricted domain of forces. An instructive way to display the supercoiled DNA 
stretching data is to plot the molecule extension versus the degree of supercoiling a (defined as the ratio of the 
number of bead turns to the number of double helix turns in the relaxed molecule), at various fixed values of the 
stretching force. Three such plots are shown on Figure |l] for typical force values. Working with molecules having a 
moderate degree of supercoiling, a < 0.1, two changes of regime occur in the force range between F = O.OlpN and 
F = 5.0pN. For forces below OAApN the extension versus supercoiling curves are symmetric under the exchange 
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a — > —<T and have been called 'hat' curves for obvious reasons. This is the domain of entropic elasticity which will be 
our main topic in the present paper. When F > 0.5 pN a plateau is developed in the underwinding region a < 0. The 
mechanical underwinding energy is transformed into chemical energy which is used to open up the hydrogen bonds 
and denaturation bubbles appear along the DNA chain. Above ApN, a plateau appears also in the overwinding region 
tJ > 0. This has been interpreted as a transition to a new structure of DNA called the " DNA-P" The double 

helix structure is clearly responsible for the ability of DNA to convert mechanical torsional energy into chemical 
energy, with a net preference for underwinding energy. 

Our work here will be restricted to the entropic elastic regime, and its aim is primarily to compute the extension 
versus supercoiling hat curves. We have used the simplest extension of the WLC model, which is still a phantom chain, 
with one single new elastic constant, the twist rigidity. We call the corresponding model a Rod Like Chain (RLC). 
We will show that the RLC model, described in terms of local Euler angle variables, is singular in the limit of a purely 
continuous chain. Its physical realization requires the introduction of a short length scale cutoff (which will turn out 
to be of the order of magnitude of the double helix pitch) . Then the RLC is able to reproduce all experimental data 
in the restricted domain of its validity which we have discussed above. The good fit to the experimental data allows 
to deduce the twist rigidity to bending rigidity ratio, which was rather poorly known so far (with uncertainties up 
to a factor or 2). A short account of some parts of this work has been discussed in a previous letter (where the 
RLC model appeared under a slightly different name: the WLRC model). We shall give here a more comprehensive 
- and hopefully comprehensible- description of our model. Our principal aim is to show how the RLC model, once 
it has been properly formulated, can be solved exactly: by exact, we mean that the predictions of the model can be 
computed by a mixture of analytic and numerical methods to any desired precision. 

The first conclusive work able to fit the experimental data from a rod like chain with two elastic constants is that 
of Marko and Vologodskii [B0|, who used a discretized model and performed some Monte Carlo simulations. On the 
analytical side, Fain et al. [25| wrote the elastic energy in terms of a local writhe formula, similarly to the one which 
we have used, but they used it only in the zero temperature limit. In an independent work which appeared just after 
ours, Moroz and Nelson ||2l[] used the same form of energy as we did, but with a different approach for handling the 
singularity of the model in the continuum limit. We shall comment extensively on the differences of these approaches, 
as well as on the obtained results. 

Developing a theory beyond the elastic regime is a much more difficult task, since such a theory should involve both 
the self repulsion of the chain (necessary to prevent a collapse of plectonemes) , and the possibility of denaturation. 
Some first attempts in these directions can be found in . 

Recently an explicit model, coupling the hydrogen-bond opening with the untwisting of the double helix, has been 
proposed It leads to a unified description of thermal and untwisting DNA denaturation in the high stretching 
force limit. 

We now indicate the organization of the paper, with some emphasis upon the "RLC model crisis" and its solution. 
In sect. H we introduce the RLC model and discuss its range of validity. We note first that the elastic energy used in 
the present paper is invariant under rotation about the molecular axis, in apparent contradiction with the double helix 
structure. We prove that added cylindrical asymmetric terms in the elastic energy are washed out by averaging upon 
the empirical length resolution Al, about three times the double helix pitch p. We also discuss the relationship to 
the Quantum Mechanics problem of a symmetric top, which appears naturally when one considers the configurations 
of the DNA chain as world lines in a quantum mechanical problem. We stress that, despite a formal analogy, there 
is an important difference in the formulation of the two problems. The RLC analog of the angular momentum is 
not quantized since, in contrast to the symmetric top, the physical states of an elastic rod are not invariant under 
rotations of 27r about its axis. This is the origin of the "RLC crisis". 

In Sect. Ill we incorporate in the RLC partition function the supercoiling constraint coming from the rotation of 
the magnetic bead. First written as a boundary condition upon Euler angles used to specify the DNA configuration, 
the supercoiling constraint is transformed, assuming some regularity conditions on the Euler angles, into an equality 
between the bead rotation angle and the sum of the twist and writhe variables of the open molecular chain. These 
variables fluctuate independently in absence of supercoiling so that the constrained RLC partition function is given 
as a convolution product of the twist and writhe partition functions. The writhe partition function Fourier transform 
is in correspondence with the Quantum Mechanics problem of a charged particle moving in the field of a magnetic 
monopolc with an unquantized charge. The spontaneous fiuctuations of the writhe provide a spectacular signature of 
the RLC model pathology: the second moment of the fiuctuation is predicted to be infinite in the continuum limit. 
This divergence comes from a singularity appearing in the RLC potential when the molecular axis is antiparallel to 
the stretching force. 

Sect. IV explains how the angular cutoff needed in order to transform the RLC model into a sensible one can be 
generated by the discretization of the chain. The introduction of a length cutoff 6, of the order A^, appears rather 
natural since an average upon has been invoked to justify our choice of the RLC elastic energy. We proceed 
next with the development of the two main tools used to get an exact solution of the discretized RLC model: a 



2 



direct transfer matrix approach and a Quantum Mechanics approach involving a regularized RLC Hamiltonian with 
a potential derived from the discretized model and now free from any singularity. 

In sect. ^ we use the two above methods in order to compute the elongatio n ver sus supercoiling characteristics of 
a long DNA chain. The results are compared with experimental data in sect. VII where the small distance cutoff b 
and the twist elastic constant are determined. 

Sect. VI presents a complementary approach, that of Monte Carlo simulations which provides the only method so 
far to introduce the constraints of self-avoidance and knots elimination. After a presentation of the results of Marko 
and Vologodskii ]3C| ], we describe our less ambitious simulations which were used to validate our analytic method. 

In sect. VIII we use the RLC model to analyse, at a force taken on the high side of our explored area, the variation 
of the torque and the writhe to twist ratio versus supercoiling. The two curves exhibit a rather sharp change of regime 
near one given value of the supercoiling angle: the torque, after a nearly linear increase, becomes almost supercoiling 
independent while the writhe to twist ratio, initialy confined to the 20% level, develops a fast linear increase. We give 
arguments suggesting that this quasi-transition is associated with the creation of plectonemc-likc configurations, able 
to absorb supercoiling at constant torque. 



Finally the sect. IX studies, within the RLC model, the extension fluctuations of an isolated DNA molecule, subject 



both to stretching and supercoiling constraints. The predictions are made for two thermodynamic ensembles: one 
at fixed torque, the other at fixed supercoiling angle. Although the two ensembles lead, as they should, to identical 
results as far as average values are concerned, they yield widely different predictions for the extension fluctuations 
above a certain supercoiling threshold (close to the one appearing in the previous section) : the fixed torque ensemble 
becomes much more noisy. We give a qualitative explanation for this peculiar behaviour and compare our predictions 
to some preliminary experimental data at fixed supercoiling. 

In sect. ^ we give a summary of the work and some perspectives on its future extensions. 



II. ELASTIC ENERGY OF THE ROD-LIKE CHAIN 



A given configuration of the rod-like chain (RLC) is specified, in the continuous limit, by the local orthonormal 
trihedron {ei(s)} = {u(s), n(s), t(s)} where s is the arc length along the molecule, t is the unit vector tangent to 
the chain, u(s) is along the basis line and n(s) = t(s) A u(s). The evolution of the trihedron {ei(s)} along the chain 
is obtained by applying a rotation TZ{s) to a reference trihedron {e^(s)} attached to a rectilinear relaxed molecule: 
e,(s)=7^(s)•eO(s). 

The rotation TZ{s) is parameterized by the usual three Euler angles 0{s), (j){s) and "ipis)- The reference trihedron 
is such that 9(s) = 0, 0(s) -I- ip{s) = loqS, where loq is the rotation per unit length of the base axis in the relaxed 
rectilinear DNA molecule. With the above deflnition, the set of s dependent Euler angles 9{s), (j){s),tjj{s) describes the 
general deformations of the DNA molecule with respect to the relaxed rectilinear configuration. It is now convenient 
to introduce the angular velocity vector r2(s) associated with the rotation TZ{s). The evolution of the trihedron {ei(s)} 
along the molecular chain is given in term of r2(s) by the set of equations: 

^e,(s) ^ (f2(s)+wot(s)) Ae,(s) (1) 

We shall use in the following the components of Cl{s) along the trihedron {ei(s)}: fli — ft{s) ■ei{s) with ei = u, 62 = n 
and 63 = t. ( The Qi are computed in terms of the Euler angles and their s derivatives in the Appendix A). The 
stretched RLC energy E^lc is written as a line integral along the chain of length L, involving a sum of three 
contributions: 

ErlC = ksT / ds {£bend{s) + £twist{s) + £stretch{s)) (2) 
JO 

The above bend, twist and stretch linear energy densities are given by: 

^ (02 sin^ e + 9^) 



(3) 

where the dot stands for the s derivative. 



£bend - ^ (^^1 + ^^2) - 



C C ■ ■ 

Stwist = —nl = —(tP + cI) COS 9f 

Fcost 



-'stretch 



-t(s).F/(fcBT) = -: 



knT 
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^bend IS proportional to the inverse of the square of the curvature radius and represents the resistance against 
bending around an axis perpendicular to the chain axis. The coefficient A is the persistence length; its typical value 
is 50nm. 

£twist is proportional to the square of the projection of the angular velocity along the molecular axis and gives 
the energy associated with a twist constraint. The twist rigidity C is known to b e in the range 50 to 100 nm; its 



determination from the experimental measurement will be the topic of section |VII . 

Sstrech is the potential energy associated with the uniform stretching force F = F.z which is applied at the free end 
of the molecule; it is written as a line integral involving t(s).z. 

The ground state configurations of the rod like chain have been studied in ||2^ . These may apply to the properties 
of small plasmids, but in the regime which we are interested in, the thermal fluctuation effects are crucial. If no 
twist dependent properties are imposed or measured, one can factor out the integral in the partition function, or 
equivalently use the C — theory, in which case one recovers the elastic energy of the Worm Like Chain model : 

- f 4 - cos^(.)F/(fc,T)) (4) 



ksT Jq \2 ds 
It is now appropriate to raise two questions about the above formulas 



A. Is a cylindrical symmetric elastic rigidity tensor adequate for the DNA chain? 

The theory which we are using here is an clastic theory which cannot be valid down to atomic scales. In particular it 
does not take into account the microscopic charges on the DNA and their Coulomb interactions in the solvent, but 
it just describes the net effect by a set of elastic constants. However, even at the level of an elastic description, one 
may wonder whether our choice of elastic tensor, which ignores the helical structure of the DNA, is valid. Indeed the 
elastic energy given by equations (|^) involves a rigidity tensor with a cylindrical symmetry around the molecular axis 
t(s). We would like to argue that such a description is a reasonable approximation if one deals with experimental data 
obtained with a finite length resolution ~ lOnm, as is the case for the experiments we shall analyse. The pitch p 
of the double helix contains approximately 10 basis which corresponds to a length of SAnm. A simple way to break 
the cylindrical symmetry is to introduce in the elastic energy linear density a term proportional to Ar2(s) = fl^ — flj- 
As it is shown in Appendix A, Ari(s) can be written in terms of the Euler angles in the following way: 

An(s) = n_L^ cos2(C(s) + + — ) 

p 

with r2_L^ = Ql + nl and C(s) = arctan ^^/(sin6'0)^ ; 

One must note that Ari(s) oscillates with a wavelength which is half the pitch p of the double helix. An av- 
erage involving a resolution length, about 6 times the oscillation wavelength, is expected to lead to a strong 
suppression of AO(s). The average Ail(s) = J dsiP{si — s)Afi(si), associated with the resolution function 
P(s) = ■^^^=^ cxp(— ifi^/A ^^), is computed explicitly in Appendix A in terms of simple Gauss integrals, assum- 
ing that the phase C(s) + and the amplitude vary linearly within the interval {s — Al,s + Al). Ignoring 
first the fl± variation and assuming that the supercoiling angle per unit length is <C wo, we arrive in Appendix A to 
the result: 

Ml{^ w An{s) exp J^-i(l!!-^)2^ 

Taking p = 3Anm, Al — lOnm one gets: ^^^^ ~ 37; it is clear the above expression of An{s) / AD,{s) is zero for all 

practical purposes. The extra term also found to be exceedingly small. The same suppression factor holds 

for rii Q2 while in the case of fii fl^ the argument of the exponential is divided by 4 but with the value quoted above 
for A^ the suppression effect is still very important. Therefore all the cylindrical asymmetric terms in the elastic 
rigidity tensor are washed out by the empirical finite length resolution averaging. 

In the same spirit there are various effects which are not taken into account by our description. An obvious one is 
the heterogeneity of the sequence: a more microscopic description should include some fluctuations of the rigidities 
along the chain, depending on the sequence of bases. However one expects that such fluctuations will be averaged out 
on some long length scales, such as those involved in the experimental situation under study. This has been confirmed 
by some more detailed study involving some simple model of disordered sequences psi. One should also notice that 
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we have not allowed the total length of the chain to vary. It is possible to add some elastic energy to bond stretching 
as well as some stretch-twist coupling in order to try to describe the behaviour at rather large forces and supercoiling 
|p^|-pT| . But as we will see such terms are irrelevant for the elastic regime which we study here. 

B. What is the precise connection of the elastic rod thermal fluctuations problem with the quantum theory 

of a symmetric top ? 

A careful inspection of the formulas (||) and (||) suggests a close analogy of the elastic linear energy density with 
the classical lagrangian of the symmetric top with A and C being, up to a constant factor, the moments of inertia. 
For given values of the Euler angles d{s), (/)(s),'0(s) at the two ends sq — and si = i of the chain, the partition 
function of the RLC model is given by the following path integral: 

Z(0i,(/.i,Vi,si|0o,'/'o,V'o,so) = Jv{9,cb,^)exp(^-^^^ (5) 

where V {9, (j), ip) stands for the integration functional measure for the set of paths joining two points of the Euler 
angles space with < < tt but no constraints imposed upon 4> and if) . Let us perform in the elastic energy integral 
over s ( Eq|^ ) an analytic continuation towards the imaginary s axis. We shall use for convenience a system of units 
where fi = c — ksT — 1 and write 3s = t. The elastic energy Erlc is transformed into —i J^^ dtC{t), where the 

Lagrangian is that of a symmetric top in a static electric field /: C{t) = ^Y^^CiVLf + f cos 9{t)^ with inertia moments 
Ci — C2 — — C and / = F/{kBT). The analytically continued partition function of the RLC model is then 

identified with the Feynmann amplitude: 

(0i,0i,i/'i,ti|0o,0o,V'o,io) - Jv{e,^,tp)exp(^iJ^ dtC{t)^ 

= (6*1, </)i,'(/'i| exp (^-i{ti ~ to)'Htop^ \0o, 0o,i/'o) 

where Ti-top is the symmetric top Hamiltonian written as a second order differential operator acting upon Euler angles 
wave functions. Returning to real value of s, ones gets the RLC model partition function as a Quantum Mechanics 
matrix element, with the substitution i(ti — to) —fSi — sq. It is now convenient to introduce the complete set of 
energy eigenstates of Htop- 'Htop \ n) ~ £'„ | n). One gets then: 

Z(6'i,</)i,V'i,si|6'o>o,V'o,so) = X! cxp-i£;„ 

n 

*„(0i,0i,Vi)*:(0o,0o,^o) (6) 

where 5'„(6', </>, ■0) = is the eigenfunction relative to the state \n). As we shall discuss later, the interest of 

the above formula lies in the fact that the above sum is dominated by the term of lowest energy. 

At this point one may get the impression that the solution of the RLC model is reduced to a relatively straightforward 
problem of Quantum Mechanics. Although this identification is valid in the case of the worm-like-chain model (WLC) 
which describes the entropic elasticity of a stretched DNA molecule with no twist constraints, it is not the case for 
the more general RLC model. The above analysis was somewhat formal in the sense that the Hamiltonian Titop was 
given as a differential operator and it is thus not completely defined, until the functional space on which it is acting 
is properly specified. Physical considerations will lead us to choose different functional spaces for the two problems 
in hands: 

In the symmetric top problem the space is that of Itx periodic functions of the Euler angles (p and -0 , but in the 
study of RLC thermal fluctuations the space is that of general functions of these angles, without any constraint of 
periodicity. 

As an illustration, let us discuss the simple case of a non-flexible rod with E = ^ksT ds'ip{s) . The associated 
quantum problem is the cylindrical rotator described by the Hamiltonian: Tirot — "^^^ -^tt periodicity of 

the motion implies a discrete spectrum for the conjugate momentum p^ = ^^"^ corresponding wave function is 

exp(zA:V'), where k is an arbitrary integer. Imposing on the molecular chain the boundary conditions ijj{so = 0) = 
, = L) ^ the partition function of the non-flexible rod is obtained by a straightforward application of formula 



Z(x,L) = ^exp (^ikx-^k'^ 
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In the experiments to be discussed in the present paper L ^ C and, as a consequence, the sum over the integer k is 
dominated by the terms fc = 0, ±1, leading to the rather unphysical resuh: 2^(x, L) — 1 + 2 cosxexp(— 2^). 

On the contrary, in order to describe the twisting of an elastic rod. The spectrum of p-^ has to be taken as continuous 
and the dicrete sum J2k ■■■ ^lust be replaced by the integral dfc/(27r)..., which yields the correct result for the 
partition function: 

Z(x,L)aexp-(^x') (7) 

Two significant physical quantities can be easily obtained from the above partition function: First, the second moment 
< X^ > relative to a situation where x is no longer constrained but allowed to fluctuate freely. Keeping in mind that 
P(x) oc Z{x,L), one gets : < x^ >= ^- Second, the torque F given by the logarithmic derivative of Z{x,L) with 
respect to x: T = ksT^. The above quantities look very reasonable from a physical point of view: 

< > scales linearly in L, as expected from the fluctuations of a linear chain with a finite correlation length C 
and the expression of F reproduces an elasticity textbook formula, if one remembers that ksTC is the usual twist 
rigidity. 

The above considerations may look somewhat simple-minded. It turns out however that they have important 
physical consequences. In order to solve the RLC model, we shall have to deal with a quantum spherical top when 
its angular momentum along the top axis is not quantized. Such a problem is mathematically singular and as a 
consequence the continuous limit of the RLC model we have considered so far will not give an adequate description of 
supercoiled DNA. We shall have to introduce a discretized version of the RLC model, involving an elementary length 
scale b about twice the double helix pitch p. This is consistent with the considerations of the previous subsection 
where the empirical length resolution A Z w 3p was invoked in order to justify the cylindrical symmetry of the tensor 
of elastic rigidities. 



III. THE PARTITION FUNCTION FOR SUPERCOILED DNA IN THE ROD-LIKE CHAIN MODEL 



A. Implementing the experimental supercoiling constraint on the pcirtition function 

The first step is to incorporate in the functional integral the supercoiling constraint which results from a rotation 
around the stretching force F of a magnetic bead biochemically glued to the free end of the DNA chain. The rotation 
TZ in terms of the Euler angles is usually written as a product of three elementary rotations, a rotation 9 about the y 
axis sandwiched between two rotations ip and (f) about the z axis, performed in that order: TZ — R{z, (p) R{y, 0) R{z, -tp), 
where i?(m, 7) stands for a rotation 7 about the unit vector rh. It is convenient, in the present context, to introduce a 
different form of TZ involving a product of two rotations written in terms of a new set of variables, 6, ijj and x = + V'- 
TZ = i?(z, x) R{v/{tp),9) where the unit vector w{'tp) lies in the x y plane and is given by: = i?(z, —tp)y. ( See 

Appendix A for a proof of the identity of the two above forms of TZ). Let us consider, before any application of an 
external rotation, the DNA segment sticking out from the bead. Assume that it is short enough so that its direction is 
not affected by thermal fluctuations. Its orientation is specified by three Euler angles, Oin,ipiniXin- If the rotation of 
the magnetic bead by n turns is performed adiabatically, the final orientation of the molecular end trihedron {ei(i)} 
is specified by the rotation: 

TZ{L) = R{z,2nn + Xm) R {vf{ipin), Oin) ■ 

We can read off easily from the above formifla: x(L) = 4>{L) + fpiL) = 27rn + Xin- It is reasonable to assume that in 
the gluing process no large surpercoiling is involved so that the angles Xin and do not exceed 2 tt. As we will see 
the relevant scaling variable for the partition function is: 

^ (8) 

In the experiments to be analyzed in the present paper, A/L ~ with 77 of the order of unity, so that Xin 
can be safely ignored in comparison to 27m. In the following we shall drop the L dependance in x(-^) and write 
X = 27rn = (piL) + ipiL). The above equation does not imply that x is a discrete variable: it is just for practical 
reasons that the experimentalists perform measurements with an integer number of turns. On the other hand, if we 
were dealing with a closed DNA chains instead of an open one, the supercoiling angle x(-^)) which is here an arbitrary 
real, would get replaced by 27r Lk where the integer Lk is the topological linking number. 
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1. A local writhe formula 

We make now an assumption which may look, at first sight, somewhat trivial but has, as we shall see, far reaching 
consequences: 

The Euler angles (j){s) and ip{s) are regular enough functions of s to allow the replacement of x by a line integral 
involving the sum of Euler angles derivatives, taken along the trajectory T of the tip of the tangent vector t(s). 

X^2TTn^cf>{L)+i:{L)^ 1^ ds (ij + <i>) (9) 

It is then convenient to introduce the total twist along the chain T^, which appears as a Gaussian variable in the 
partition function. 

ds ~ dsila — 
Jo Jo 



r„ = / ds^ dsn^^ / ds + cos (lo) 



where i^sis) is the projection of the instantaneous rotation vector $7 onto the unit tangent vector t(s), given in terms 
of Euler angles in Appendix A. We now define a 'local writhe' contribution xw by substracting the total twist 
from the supercoiling angle x 

Xw=X-Tu<^l </>(!- COS0) (11) 







The above decomposition follows from our analysis of the empirical supercoiling angle and a mere line integral 
manipulation. This is reminiscent of the decomposition of the linking number into twist and writhe for closed chains 

HQ- 

Formula (|ll) has been first obtained by Fain et al. by adapting to the case of open strings a formula first 
derived by Ful er for closed strings. We have decided to include an explicit derivation here in order to be able to 
describe some subtleties in its use, which will require some regularization procedure. Since t(s) is a unitary vector, the 
curve r lies upon the sphere 5*2 and it is parametrized by the spherical coordinates (6'(s), (/)(s)). This representation 
is well known to be singular at the two poles 9 = 0,9 = tt ( the choice of the z axis is not arbritrary but dictated 
by the experimental conditions: it coincides with the rotation axis of the magnetic bead and gives also the stretching 
force direction). With the usual convention, 9 is required to lie within the interval [0,7r]. This condition plays an 
essential role in the quantum- like treatment of the WLC and RLC continuous models, to be discussed later on. When 
the trajectory F passes through the poles the restriction imposed upon 9 implies that the function (j){s) suffers form 
a discontinuity of ±7r and this invalidates our derivation of formula (pi]). To get around this difficulty, we have to 
pierce the sphere 5*2 at the two poles. The two holes arc defined by two horizontal circles with an arbitrary small 
but finite radius e with their centers lying on the z axis. Note that 5*2 has now the topology of a cylinder so that 
the piercing of the sphere ^2 provides the necessary topological discrimination between the two distinct physical 
states {9{s), (j){s),^{s)) and (6'(s),0(s) + 27rn, V'(s)). A careful evaluation of formula (pT|), involving continuously 
deformed F trajectories drawn on the pierced sphere 52 in order to avoid the poles, gives results in agreement with 
those given by the non local closed loop writhe formula j2^, while a naive application of (|l^) would lead astray. 

2. Solenoid and plectoneme supercoil configurations 
To illustrate this point, we shall evaluate xw for right handed solenoid and plectoneme supercoil configurations, 
having arbitrary orientations. Let us begin with the simple case where the solenoid superhelix is winding around the 
z axis. The tangent vector tip trajectory Tgoi consists of n turns along the horizontal circle cut upon the unit sphere 
^2 by the horizontal plane z = cos 6*0 ( tan — where R and P are respectively the radius and the pitch of the 
superhelix.) The writhe supercoiling angle Xw{^soi) is easily seen to be given by 27rn (1 — cos^o)- 

The plectoneme supercoil configuration is composed of two interwound superhelices connected by a handle. The 
helices have opposite axis but the same helicity. The plectoneme t trajectory Fpjec is written as the union of three 
trajectories: Fpjec = ^\oi ^ ^han U F^^;. F^^; is just the right handed solenoid considered above. Than is associated 
with the handle connecting the two superhelices and its contribution to xw can be neglected in the large n limit. F^^^ 
is the tangent vector trajectory generated by the superhelix winding down the z axis; it is given by performing upon 
the F^^j the transformation defined by the change of spherical coordinates : 4>' — —(p , 9' ~ tt — 9. One readily gets 
XwiTpiec) by writing : xwi^piec) = Xiv(rL;) + XwiTl^i) = -47rn cos 6*0. 

We proceed next to the deformation of t trajectories xwi^soi) and Tpi^c by applying a rotation about the y axis of 
an angle a which is going to vary continuously from to tt. We will find that a proper use of formula (^l]) leads to 
rotation invariant results in the limit of small but finite angular cutoff. As a first step, let us show that xw{^) can 
be written as the circulation along F of a magnetic monopole vector potential. This technical detour will also be used 
in the forthcoming derivation of the RLC model Hamiltonian. 
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3. The local writhe xw o,s a magnetic monopole potential vector line integral 
Our aim is to prove that the writhe supercoiled angle xw can be written as the Une integral J d(j) = f Am (r) dr 
where = (1 — cos 9) will be identified as the (j) spherical component of the potential vector Aj„(r) of a magnetic 
monopole of charge unity. Following the well known Dirac procedure, we write Am(r) as the potential vector of a 
thin solenoid of arbitrary long length, the so called Dirac String £, lying along the negative half z-axis : 

A,„(r) = / dl A V,^- = ^^^^ (12) 

where I = uz with ~oo < u < 0. Using the above expression of A,„(r), one derives the two basic equations which 
allow the writing of xw the circulation of the magnetic monopole potential vector, valid when r + z • r > : 

A„(r) • dr = (1 - cos6i)d(/) (13) 
B™(r) ^ V A A,„(r) = ^ (14) 

where Bm(r) is indeed the magnetic field produced by the magnetic charge unit. 

Let us now define S+ (F) ( resp S- (F) ) as the part of the sphere S2 bounded by the closed circuit F which is run 
anticlockwise ( resp clockwise ) around the surface normal. ( Note that 5'+(F) U 5'_(F) — S2 )■ Let us call Se the 
part of iS'2 defined by the south pole hole and assume that 5'+(F) n Eg = ( as a consequence: 5'-(F) n Se = Se). 
Applying the Stokes theorem to formula ( |lj ), we can write: 



Xwir) ^ f Am{r) ■ dr ^ dS • B„,(r) = ^(5+(F)) (15) 

J r J J s+(r) 

where yl(S'+(F)) is the area of the spherical cap S'+(F). The closed circuit F is assumed here to be 
run only once; if it is run n times, as in solenoid and plectoneme configurations, ^(iS+(— )) should 
then be multiplied by n and the formulas obtained before are immediately recovered. A similar for- 
mula, valid mod 2tt, has been derived previously by Fuller j2^, in the case of closed molecular chain. 

4-. On the rotation invariance of the writhe supercoiling angle xw 
We are going now to follow the variation of xiv(F) when the initial tangent vector t trajectory Fq is moved contin- 
uously upon the sphere by applying a rotation about the y axis by an angle a, varying continuously from to tt: 
Fq F(q!) = R{y, Oi)To. We shall choose Fq = Tgoi but the following analysis can be easily extended to plectonemes 
or more general configurations involving closed t trajectories. 

The crossing of the north pole is harmless because 1 — cos 9 vanishes at 6* = and everything goes smoothly until 
F(a) meets the south pole hole { 9 = n ) when a is approaching tt — ^o- When a < tt — the south pole stays outside 
S'+(r(a)) and xw{r) is given by n times the spherical cap area: A{S+{r{a)) = ^^(^^.(Fo)) = 27r(l — cos 6*0). In other 
words xw(r(a)) = XwO^o) ii a < n - 9q. 

In order to avoid the hole F(a) must be continuously deformed for a > n ~ 9o into the circuit T'(a) — F(a) U Fg 
where the path F^ consists in n turns run anticlockwise around the south pole hole of radius e. (We ignore the 
two ways narrow lane connecting the two loops which gives a vanishing contribution.) It is clear now that S'-|_(F(q!) 
contains the south pole, so we must apply the Stokes theorem to the clockwise spherical bowl S'_(F(q;)): 



XwiT{a)) = -n / dS ■ Bm{r) = ~nA{S-(T{a)) = n{A{S+(T{a)) ~ 4tt) 

J J S-{r{a)) 

The contribution of the path F^ is readily found to be 47rn(l + 0(e^)) so it just cancels the —Airn term in the 
above result. As a consequence, the writhe supercoiling angle associated with the south pole avoiding t trajectory 
Xw(r'(a)) coincides when n — 9() < a < tt with the a = initial result: xw(X'i<^)) = Xwi^a) — 27rn(l — cos 6*0) up 
to corrections O(e^). 

In conclusion, we have shown that a proper use of formula ^iJj), in the simple case of solenoidal configurations, 
leads to rotation invariant xw *^ the limit of small e. 

This result can be easily extended to the plectonemic configurations and even to more general ones. The crucial 
cancellation is taking place in the vicinity of the south pole crossing and thus does not depend on the detailed shape 
closed trajectory Fq in the t space. If one adds the minimal extra string section to the solenoid or plectoneme, 
necessary to get a closed loop in the physical space, our findings are in agreement, in the limit of large n, with those 
obtained from the non-local closed loop writhe formula pSf, which is explicitly rotation invariant from the start. 



Formula (11) is thus a correct description of the problem provided the trajectory is obtained from the straight line 



{9[s) — 0, 4>{s) = 0) by a continuous transformation on the unit sphere pierced a,t 9 — n. A similar condition was 
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obtained by Fuller |24| in his local Writhe formula for closed strings. An alternative derivation of our result, more 
similar to the approach of Fain |25|, could be to close the open string by some straight line closing at infinity and rely 
on Fuller's formula, but this procedure raises subtle questions about the contribution of added pieces which made us 
prefer the present explicit derivation. 



B. The Rod Like Chain Hamiltonian 



The partition function for a fixed value of x is given by the path integral in the space of Euler angles : 

Z(x, F)^ j Vie, ^,^)s(^x- ds{^ + 7^)^ exp (16) 

The functional space 6{s), (f){s),ip{s) should be specified in accordance with the previous considerations, in particular 
the tangent vector t paths must bypass the holes on the pierced sphere 52. This can be achieved by introducing 
explicitly in Eulc a short range repulsive potential near 9 = n. 

Two extra constraints, non local in the tangent vector t space should have been, in principle, implemented in the 
above partition function. 

The first one concerns self avoidance effects induced by the screened Coulomb repulsion between different molecular 
sections. Even the simple version of this constraint, used in supercoiled plasmids Monte Carlo simulations is 
exceedingly difficult to implement in the Quantum Mechanics inspired formalism of the present paper. Here Coulomb 
induced self avoidance wail be ignored as all quantitative analytic approaches did so far, including the celebrated 
WLC computation. It will appear that these effects play a limited role in the low supercoiling regime ( \a\ < 0.02) 
to be analyzed in the present paper, a regime and experimental conditions where the WLC model is also working 
beautifully. 

The second non-local constraint not implemented in the present work concerns the exclusion of knotted chain 
configurations. In principle, knotting an open DNA chain is not forbidden, but such a transition is expected to be 
inhibited on the time scale of the actual experiment since the contour length of the molecule L is only 1.7 times the 
circumference of the bead glued at the free end of the molecule. 

Using a Fourier representation of the Dirac S function we can write the partition function as a Fourier integral: 
Z{x,F) — J dk exp(—ik x)Z(k) where the Fourier transform Z(k) is given by the functional integrals: 

Zik) = J Vie,^) exp (^^kxw - ^^£!^1±|^^ Z^k) 

Zrik) = j Vii,) exp [ikT^- (17) 

where the above factorization follows from the identity: ds{(j) + V') = ^tu + Xw, using the elastic energy densities 
defined in (^). The Gaussian path integral upon is readily performed and gives the Fourier transform of the twist 
partition function of equation (^ Zrik) = exp(— ^^), up to a trivial constant. The partition function Fourier 

transform Z{k) can then be written as the product: Z{k) = ZT{k)Zw{k) where Zw{k), which is interpreted as the 
writhe partition function Fourier transform, is given by a path integral upon 9 and 0, with the effective energy: 

Ew{k)^^^ + ikxw (18) 

Using the results of the previous section ones sees immediately that i k xw can be written as i times the line integral 
J Am(r) dr where A„j(r) is the vector potential produced by the magnetic monopole, having now a charge k. If one 
performs, as in Section II. B, an analytic continuation of the integral over s in Ew{k) towards the imaginary axis the 
i factor disappears in the writhe line integral. One arrives to the action integral of a unit charge particle moving on 
the unit sphere under the joint action of an electric field / and a magnetic monopole of charge k. In order to get the 
corresponding Hamiltonian Hjn^cik) we apply a standard Quantum Mechanics rule used to implement the switching 
on of a magnetic field: we replace in the kinetic terms ij^P^, the particle momentum p by p — eA„i{r) (here e=l). 
The term linear in Am(r) disappears by averaging over the final angle = 4>iL) and we are left with the diamagnetic 
term: 

1 2 / \ k"^ I — cos9 
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Adding to the WLC Hamiltonian, one arrives to the dimensionless Hamiltonian HRLc{k) associated with the 
partition function Zw{k): 

1 d . „d „ k'^ I — cos 9 

2 sm & o& oO 2 1 + cos 

where we have taken, as the length unit, the persistence length A and introduced the dimensionless force parameter 

" ksT- 

Let us stress that the above Hamiltonian is a straightforward consequence of the expression for the elastic energy ( 
Eq and our expression ^ giving the supercoiling constraint x in terms of Euler angles derivatives. 

Moreover, as it is shown in Appendix B, HRLc{k) can be derived directly, without having to introduce the local 
writhe given by formula The method involves the quantization of the symmetric top Hamiltonian as a differential 
operator acting upon non-periodic functions of the Euler angles 4> ^^nd tp. This is required to describe the thermal 
fluctuations of an clastic rod ( see section II. B) and k appears as the continuous angular momentum along the top 
axis. 

Introducing the eigenstates and the eigenvalues of Hrlc, HRLc'^n{k,9) = e„(Q;, /c^)^'„(A:, 0), the Fourier trans- 
formed partition function Z can be written as the sum; 

Z = ^vl/„(fc,0o)*„(fc,^L)exp('-^ U,{a,e)+^-^\\ (21) 

(Note that Hrlc being a real operator the eigenfunctions can always be taken as real.) 

In the large L limit, the sum over the eigenstates is dominated by the one with lowest energy eo(a, A;^), if L/A ^ Ae 
where Ae is the energy gap between the ground state and the nearest excited state of Hrlc- This gives the approximate 
expression for the partition function Z{x,F) : 



Z~ J dk cxp(^-j(^eoia,k^) + ^^ -ikx^ (22) 



To get the above equation, we have also neglected in Z{k) the prefactor involving the ground state wave function: 
^'o(A:, 6'o)'I'o(fc, 6*^). It leads to finite size corrections of order A/L, which are in general dominant with respect to 
those associated with the excited states, which scale as exp{—AeL/A). 

C. The pathology of the RLC model in the continuous limit. 

In order to comply with the prescription given for the path integral ( p^ we should have added to HRLcik) a strong 
repulsive potential acting in the interval: tt — e < 6* < tt, in order to avoid the 6 = tt singular point. We shall refer to 
the e — + limit as the 'continuous limit', and show that in that limit the model is pathological. 

Let us derive from the ground state energy eo(a, /c^) of the Hamiltonian Hrlc, some observable properties of a 
long RLC. In the following the variable z = k"^ \s assumed to take positive real values. 

If instead of constraining x measures its thermal fluctuations, their probability distribution is just P{x) Z . 
Let us consider the second moment as we did in subsection II. B for the case of a non-flexible rod. 

2 1 d^Zik) 

C A fc2^o ^ ' 

We recognize in the first term of equation ( |24| ) the contribution to < > from the twist fluctuations, ^, obtained 
previously. The second piece of ( p^ gives the contribution from (xvk)' turns out to be divergent. Evaluating 

eo{a, fc^) at small fc^ from standard perturbation theory, we flnd (xw^) = (^M)(*J'ol(l ^ cos 6*) /(I + cos0)|4>o) where 
^o{9) is the groundstate eigenfunction of the WLC Hamiltonian (which is Hrlc at fc = 0). As ^o{6 = tt) (for any 
flnite force), we get clearly a logarithmically divergent result. To study analytically this peculiar theory two methods 
have been followed, associated with different limits. 
i ) The weak force limit: 

When a = the eigenvalue problem for Hrlc can be solved exactly. The eigenfunctions are given in terms of 
Jacobi polynomials : 

*„(fc,6i) = (l + cos6i)'=P^"^2fe'(cos6i) 
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and the associated eigenvalues are: 

eo{a = 0, fc^) = {{2n + l)k + + n) /2 
In the presence of a small force (a ^ 1), the ground state energy shift is easily calculated to second order in a: 

Note that eo(a, fc^) is not an analytic function of z = fc^. 

iijThe small k limit: In order to build a perturbative expansion near = it is convenient to factorize in the wave 
function the singular term near = tt by writing : ^Q{k,9) = (1 + cos0)'^xo(fcj The wave function Xo{k,9) obeys a 
new wave equation which is amenable to a well defined perturbation expansion in k. It will turn out that the values 
of k relevant for the evaluation of Fourier transform are of the order of A/L, so that a first order perturbation is 
sufficient (at least when a is small enough in order to guarantee that ^q{7t)^ stays of the order of unity). We have 
found that the ground state energy is given to first order in k by: £0(0;, k"^) — ^0(0^, 0) + k^o{TT)^ + 0{k^). It can be 
checked that the perturbation result of method i) is consistent with the above equation. The Fourier integral can be 
readily performed, leading to a Cauchy type probability distribution, which is best expressed in terms of the reduced 
variable 77 = xA/L: P{rj) ~ i <I>o(7r)^/(<I>o(7r)'* + 77^) and thus leading to a diverging second moment. This last point 

can also be seen from the fact that ^^"g^Ji ^ °^ i the small k limit. The relative extension of the chain in the 
direction of the force is given by {z) / L = {A/L) ^g^^ . The partition function in the present case is easily found to be 
: Z{'q,a) oc exp eo(a, 0)) P{ri). It leads to: {z)/L = -^(a,0) + 0{A/L), i.e. the same resuh as the WLC 

model. 

The present computations show that the continuous RLC model leads to predictions which are both pathogical ( 
infinite second moment of the writhe spontaneous fluctuations ) and in .striking contradiction with experiment ( absence 
of variation of the average extension with respect to the supercoiling angle x moderate forces ). All these peculiar 
features have been confirmed with explicit computations using a regularized Hrlc hamiltonian. We have indeed found 
that < Xw ^ scales like ^ log ^ and that the extension versus supercoiling curves flatten up in the small e limit. 

As a last remark we would like to note that if the variable z is continued analytically towards negative value 
z —k} the groundstate energy develops an imaginary part. This a clear indication that -ffflLc(*'^) has no stable 
ground state. 

IV. DISCRETIZATION AS A REGULARIZATION METHOD FOR THE ROD LIKE CHAIN MODEL 



The introduction by hand of an angular cutoff near 9 = it appears as an ad hoc procedure. We are going to show 
in the present section that a discretization of the chain involving an elementary link of length b generates an angular 
cutoff: sin^ 9 > We have seen in section II. A that the RLC model with its cylindrical symmetric rigidities tensor is 
realistic only in presence of a finite resolution A I in the length measurement. The existence of a length cut-off 6 ~ A Z 
appears not only natural but necessary in the present context. The non trivial fact is that the "mesoscopic " elastic 
properties, taking place on the length scale of the whole supercoiled molecule (about ten microns), are sensitive to 
the existence or not of a cutoff in the range of few nanometers. We should also stress that length cutoff effects are 
found to be reduced to the level of few percents when supercoiling is absent. Therefore in that case the continuous 
version of the WLC model remains a remarkably good description of the DNA force extension measurements. 

Moroz and Nelson |2^] have devised a computation procedure where no angular cutoff is introduced from the 
start. They consider the situation where the torque acting upon the molecular end is kept fixed so that the relevant 
Hamiltonian is just if hlc (*'*)• The pathology discussed above is expressed in the fact that this Hamiltonian has no 
stable ground state. In order to deal with such a situation, they use a perturbation method near the high force limit 
where the relative elongation is close to unity. More precisely, they construct a perturbation expansion involving 

negative power of K =^ y^a — ^k^. The lowest order Hamiltonian potential is given, up to a constant, by the small 

angle approximation: ^K'^9'^ = i(a — iK^)0^. Since Hrlc{'^k) is not a self adjoint operator, the series diverges and 
it is expected to be, at best, an asymptotic series. The regime of validity of the prediction is difficult to assess in this 
approach since one is basically approximating a divergent ground state energy from the first terms of the perturbation 
series. As long as one sticks to a finite order, the 9 = tt singularity does not show up. In order to make contact with 
experiment, Moroz and Nelson have been forced to restrict their analysis to the domain of force and supercoiling where 
K > Kc = 3. In our approach, we have instead regularized the model explicitely, as we shall discuss. This allows to 
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get a better control of the theory and it turns out that our results hold on a much wider range of force-supercoiling. 
It can even be compared in the F = limit to supercoiled plasmids experimental data and Monte Carlo simulations. 



See Section VIIl 



A. Construction of the transfer matrix associated with the discretized RLC model 

The continuous elastic RLC is transformed into a chain composed of N elementary links of length b. The continuous 
arc length s is replaced by a discrete set: s„ = n6 with < n < iV and L = Nb. We have now to specify the 
discretization rules. For the bending linear energy density Sbend, there is a rather natural prescription: it is to replace 

\^\'hy±\K-K-if. 
We are led in this way to the first discretization rule: 

£bend =^ ^[1 - C0S(6'„ - 9n-l) + (1 - COs((/)„ - SmOn Sm0n~i] (25) 

The above discretized linear density is endowed with remarkable feature: it is 27r-pcriodic with respect to the angular 
finite difference — (t>n-i- This property plays an essential role in the derivation of the transfer matrix T{9„, 9n-i) 
which allows a direct computation of the partition function Z{k), averaged upon the azimuthal angle 4'{L). 

For the term generated by the supercoiling constraint: ik{l — cos9{s))(j), we shall impose the same periodicity 
condition, together with a symmetry under the exchange 0„+i <^ 0„, in order to insure the hermiticity of the transfer 
matrix. 

We arrive in this way to the second discretization rule: 

a\ 1 V 1 • J. \(^ cos0„ + cos6'„_i\ 

(1 - cos^j (j) =^ -sm((/>„ - 0„_i) II 1 (26) 

The discretized version ZN{k, 0jv, c/jn) of the writhe partition function can be written as the following 2N dimension 
integral: 

ZN{k,9N,(l>N) ^ / Y\_^4'n-ld{cOs9„-i) exp-b {A£{k,9n,9n-l,(l)n - (t)n-l)) 



where A£{k, 0„, 0„_i, </>„ — (t>n-i) is the discretized linear energy density obtained from (18) by the replacement rules 
( p5| ) and (p6|). It is convenient to introduce the partition function Zn{k,9n,(j>n) corresponding to the intermediate 
value Sn — bn. One gets then the recurrence relation: 

Zn{k,9n,<j)n) ^ J d(l)n^l d {cOS 9n-l) CXp -6 {AS {k , 9n , 9^-1 , (/in - (t^n-l)) 

Z„_i(/c, 6'„„i, 0„_i) (27) 

Since we are ultimately interested in the partition function Z{k) averaged upon the azimuthal angle (t>{L), we shall 
define : 

2„(fc,0„) — - — / d(f)n Zn{k,9m (j^n ) (28) 

where M is an integer which can be arbitrarily large. We perform the above integral upon 0„ in the two sides 
of equation (p7|). In the right hand side we change the integration variables 4'n,4>n-i into the new set: u„ = 
4'n — 4'n-i,4'n-i- Using the periodicity of the integrand with respect to m„, we can make the replacement: 

/ du„ . . . =^ — / dur, ■ ■ ■ 

2TrM J_^,j_^^_^ 2nJ_, 

The integral upon ii„ can then be performed exactly in terms of the Bessel function of second type Io{z) and we arrive 
finally to an explicit recurrence relation involving the i/i-average partition function Zn{k, 9n) : 



Znik,9n)= Sm9n-ld9n-lT{9n,9n^l) Zn-l{k,9n^i) (29) 

Jo 

T{9i,92,k^) ^exp~{j (1 - cos(6'i - 6*2) + sin^i sin6'2) 

+ 1^ (cos 01 + cos 02 ) } /o (/ (01 , 02 , fc' ) ) (30) 
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where the Bessel function argument is given by : 



To get some basic elements of the transfer matrix formalism it is useful to rewrite the recursion relation ( p9[ ) within 
an operator formaUsm: 

\zn) = T{a, where T(a, fc^) is the operator associated with the transfer matrix T{9i, 6*2, k^)- Let us define 

as \ti{a,k'^)) the eigenstate of T{a,k^) associated with the eigenvalue ti{a,k'^). Performing an expansion upon the 
above set of eigenstates, we can write \zn) as follows: 

\zn) = r(a, \zo) = J2 ^')'^ k^)){U{a, k") \zo) 

i 

Let us assume that N is large enough so that the above sum is dominated by the contribution from the lowest 
eigenvalue. The partition function Z{x) can then be written as: 

= j dk exp^-^+ifcx^ io(a,A:^)^ 
It is now convenient to introduce the new definition : 

eo{a,k^)^ -^lnto{a,k'') (32) 



Note that we had already defined eo{a, fc^) as the lowest eigenvalue oi Hulc', it is easily verified that the two definitions 
coincide in the limit b/A <^ 1 modulo an irrelevant constant. With this new definition the equation ( |2^ ) giving Z[x) 
holds true for the discretized RLC model, within the transfer matrix formalism. 



B. The angular cutoff induced by the discretization of the RLC model 

In this section we are going to discuss what turned out to be a crucial milestone in our work: the understanding 
that a discretization of the DNA chain could provide the angular cutoff needed to make sense to the RLC model. The 



first indication came, in fact, from MonteCarlo simulations which will be discussed later on in section VI. We would 
like to show here that, indeed, the angular cutoff comes out from the transfer matrix solution of the discretized RLC 
model. Moreover our analysis will lead us to the formulation of a regularized version of the continuous RLC model. 
Our procedure involves an analytic computation of the partial derivative with respect to k^ of the transfer operator 

ground state energy eo(a, /c^), as defined by equation ( p2[ ) : dk^eo = ^'^°q^2 ■ This derivative is very important phys- 
ically, for the following two reasons. First, we remind that limj,2_,o 9j,2eo gives the second moment of the spontaneous 
"writhe" fiuctuations {Xw^) (see equation ( p^ ) ), which was found to be infinite within the continuous RLC model. 
Second, for imaginary finite values of fc = in, dk^eo enters in an essential way in the parametric representation of the 
extension versus supercoiling curves, to be derived in the next sections. 
The partial derivative 9^2 cq is expressed, using Eq. (|3|) , as 

deo{a,k^) _ ~A dto{a,k^) _ -A df{a,k^) 

bto{a,k^) (9fc2 6to(a,fc2) dk^ ^ °' ^ ' 

The par tial derivative with respect to fc^ of the transfer matrix T{0i,d2,k^) is easily obtained from equations (jsT 
and 



dT{9i,e2 ,k 



T{0i,02,k'^)W{0i,92) 
1„ ...(1 



V 2 



cos 62 

WiOM ^ -',W0.e2,k^} r ^^^J^^^,^' (34) 



We have defined the small x cutoff function 
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Rix) = ^ (35) 

which behaves Hke x when x ^ 1 and goes rapidly to 1 when x > 1. The next step is to investigate the variation of 
T{9i, 6*2, fc^) as a function of A9 — 62 — 61 for a fixed value of 6 — {6i+62)/2 when A/b^ 1. For the sake of simplicity 

we limit ourselves to value of such that ^ 1 so that f {61,62, k"^) can be approximated by ^ sin 6*1 sin ^2- 

The transfer matrix is then given by the somewhat simplified expression: 

T{6i,62,k^) = exp-{4 (1 - cos(0i - ^2)) + ^ (cos 0i + cos 02)}So(4 sin 0i sin02) (36) 
2 A b 

where we have introduced the auxiliary function: Bo{x) = exjp{—x)Io{x) which is slowly varying with x: constant 
near the origin it behaves like l/^/x for large x. This formula shows that T{6i,d2,k'^) is strongly peaked at small 
values of |A6'|, of order y/b/A. This suggests the decomposition of W{6i, 62) as a sum of two terms involving a local 
potential Uw{6) plus a small correction: 



W{6^,62) = {UM + Uw{62))/2 + AW{6u62) (37) 
UU0)-Wi6,6)^-±Ri^si.^6){l^) (38) 

A Taylor expansion of AW(6'i, 6*2) with respect to A6, with 6 kept fixed, gives: AW oc A0^ (l + 0(A6'^)). It confirms 
that the contribution of AW relative to that of W is of the order of b/A. To be more quantitative, we have computed 
the two averages (W) and (AW) using as weight function T(0i, ^2) sin^i sin02- Taking b = AAA we have obtained 
|(AW)/(W)| = 0.013. It appears then quite legitimate to set AW = 0. This approximation allows us to write the 
partial derivative of the transfer operator T as a symmetrized operator product: 



af(a,fc2) {f{a,k^)U^ + U^f{a,k^)'^ 



dB 2 
where Uw stands for the operator associated with the local potential Uw{6). 



(39) 



Introducing this formula in eq. (33) and remembering that \to) is an eigenstate of T, the partial derivative eo{a, k"^) 
takes the simple form: 

aeo(a,fc^) _ -Ait,\fu^W\to)) _ 1^^^, n_<^\ ^^^,^.e)\t,) (40) 



2tQ{a,B) 2^"\l + cos6J H 

The above result provides a very clear evidence that the discretization of the RLC model generates the angular 
cutoff around 6 = n which we needed. Indeed, the integral involved in the above quantum average is now well defined 
for 6 — 71, due to the presence of the cutoff function R{j- sin^ 6). 

Formula ( ^ ) leads in a natural way, to a formulation of a regularized version of the continuous RLC model. The 
regularized operator Hamiltonian -ffj^xc obtained from the Hamiltonian Hrlc given in equation ( pO| ) by multiplying 
the singular " writhe " potential by the cutoff function R{j- sin^ 6). Standard Quantum Mechanics rules applied to 
^RLC S^^^ 9f^2eQ a formula identical to (^), provided it is legitimate to neglect the small difference between jig) 
and 1 5*0), respectively eigenstate of T and H^j^^u{k). We have verified that it is indeed so by comparing the results 
obtained from the above regularized RLC continuous model and those given by the transfer matrix method: they do 
coincide to the level of few %, at least in the domain of stretching force explored in the present paper. 
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FIG. 1. The elongation versus the reduced supercoihng parameter o — jp— where L^o stands for the number of double helix 
turns in absence of external constraints. The curve F — 0.2pN, a typical "hat" curve, corresponds to the regime of entropic 
elasticity which is well described by the RLC model introduced in the present paper. The curve F = IpN exhibits a plateau in 
the underwound region (cr < 0) which is associated with the denaturation of the DNA. In the third curve F — 8pN a second 
plateau has developed for (a > 0); it has been interpreted as an induced transition to a new structure of DNA: the "P-DNA". 



V. THE RLC MODEL FOR A QUANTITATIVE ANALYSIS OF SUPERCOILED DNA STRETCHING 

EXPERIMENTS 



In this section we give the basic tools allowing to compute, within the RLC model, the various quantities measured 
experimentally. 

We shall first study the "hat" curve which is a graph giving, for a given force F, the relative extension of molecule 
(^^) ^-long the stretching force F, as a function of the supercoiling angle x = 27rn. Figure |l| gives three examples 
of such "hat" curves taken within a rather large range of forces. The supercoiling angle is parameterized here by the 
ratio fj = where L^o stands for the number of double helix turns in absence of external constraints. Positive values 
of a correspond to overwinding, negative values correspond to underwinding. As it is apparent in Figure |l|, only the 
curve associated with a force of 0.2 pA^ looks really like a "hat". That is so for all the extension versus supercoiling 
curves taken in the force range .06 pA^ < F < 0.45 piV. In contrast, when F > .5pN a plateau develops in the 
underwinding region. This suggests that the torsional energy is converted into chemical energy instead of entropic 
elasticity energy. More precisely, convincing experimental arguments have been given in favour of the following 
mechanism: the underwinding energy is used by the molecular system to open the hydrogen bonds linking the bases 
and, as a consequence, denaturation bubbles appear along the DNA molecule. When the force is further increased, 
say above 3piV, a plateau appears also in the overwinding region. This has received a very interesting interpretation 
PI as a transition towards a new structure of the DNA: the so called " DNA-P " , which is also predicted by numerical 
simulations. 

In the present paper, we shall focus our analysis on the true " hat " curves, i.e those symmetric under the exchange 
a <^ — (7, as observed for F < 0.45 pN. 
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A. The saddle point method 



The experiments suggest that the variations of the relative elongation (^^) versus the supercoiling angle x scale 
as a function of x/L- It is convenient to introduce an intensive supercoiling variable rj = X-^/L which will turn out to 
be at most of the order of few units in the domain we are going to explore. It is related to the usual variable a by : 

r?=^ = — a = 94.8a (41) 
L p 

where we have used the values for the pitch: p = SAnm and the persistence length: A — bl.Snm. 
Let us rewrite the partition function Z(x, F) given by equation (^2|) in terms of scaling variables: 

Z{r^,a)~ / dfc exp-- f eo(fc^a) + — -iry/cj (42) 

The above partition function can be computed by the saddle point method in the limit L/A ^ 1 with rj kept fixed. 
The saddle point is imaginary, fee — iK{a) and is given by the equation: 

C+2^("'--)-- (43) 
The saddle point contribution to the partition function ^(77, a) reads as follows : 

ln{Z{rj,a))^-j (^eoia, ^k') - ^ + v + 0{1) (44) 

Let us first compute the torque F acting upon the free end of the molecule. The experiments to be analyse in this 
paper were not designed to measure F but there is an experimental project at ENS-Paris aiming at its direct empirical 
determination. 

F dlnZ AdlnZ 



ksT dx L di] 



dn f kA ^ deo ^ 2 



= .+ ^^,- — -2.^(a,-.^)J=. (45) 

where the term proportional to ^ vanishes because of the saddle point equation (|4^). Therefore we have found 
that K = 3( fcc) is equal to the torque F in units of ksT. One can introduce the thermodynamic potential — 
— In {Z{ri, u)) — Kx and verify that the supercoiling angle % given by equation ( ^ ) satisfies the thermodynamic relation: 
y = f-S-^l 

Ok \kBT J ■ 

We are now going to compute in the same way the relative molecule elongation {^^^) : 

{z{L)) _kBT d\n Z _ def){a,-K^) 



L dF da 



(46) 



As in the computation of the torque the term proportional to ^ does not appear because it is multiplied by a 
factor which vanishes at the saddle point. In other words, the elongation is given by the same expression whether 
the experiment is performed under the conditions of fixed surpecoiling an^le x or fixed torque F. In contrast, the 
situation is very different for the elongation fluctuations: ( z{LY ) — { z{L) ) turns out to be much larger if measured 
when the torque is kept fixed instead of the supercoiling angle, as wc shall discuss later. 

From the knowledge of eo{a,—K'^), using jointly equations (^Jj and one gets a parametric representation of the 
" hat " curves, the parameter being the torque in ks T unit. 

B. Solving the quantum mechanics eigenvalues problems associated with the RLC model 

In order to complete the description of the theory we now sketch the methods used to get the final theoretical 
ingredients: the groundstate eigenvalue eo{a, fc^) and its partial derivatives. We have followed two approaches. 



16 



1. Method a): Ground state eigenvalue of the hamiltonian H^ij^c associated with the regularized continuous RLC model 



In the first approach, the computations are performed within the continuous RLC model regularized according to 
the cutoff prescription derived in section IV. B. What we have to do is to solve the ground state eigenvalue problem: 
^RLC = eo^'o(^)- This ground state wave function is obtained as an appropriate solution of the ordinary 

differential equation associated with the eigenvalue problem: 

{sine ^oie)) + {~Vr{a,-^^) + eo) ^oiO) = 0. (47) 



2 sine 80' 8e 
where the regularized potential Vr{a, k'^) is given by 

,2n . fc2 1-C0S6I /i(4sin2 6i) 

Vria,k^) ^ -acos0 + —- -^^^ ^ 48) 

^ ^ 2 1 + C0S6I lo{^sin^0) ^ ' 

We search for a solution which satisfies regularity conditions both for — and = n. These requirements are 
necessary and sufficient to guarantee that the differential operator H^j^f-, is a self-adjoint operator. They can be 
fulfilled if and only if the reduced energy e belongs to a discrete set of values. As a first step, one constructs by 
series expansion two solutions of (^), '^a{0) and ^6(0), which are respectively regular for — and = tt. One 
then proceeds to an outward numerical integration of ^'a(^) and inward numerical integration of ^6(^) up to an 
intermediate value of = Oq. The energy eigenvalue equation is obtained as a matching condition for the two wave 
functions, which ensures the regularity of the eigenfunction for the whole physical domain of 0. One writes, for = 0o, 
the equality of the logarithmic derivatives of '^a{0) and '^b{0) ■ 

Equation (^) is then solved by a standard iteration method requiring a trial approximate eigenvalue. The construction 
of the parametric " hat " curve for a given value of a begins with a small value of k ( the energy £9(0;, 0) is easily 
obtained from the results of ref . ) . The trial eigenvalue needed to solve eigenvalue equation ( ^7|) for a small value 
of K is easily obtained from a first order perturbation calculation. Since the present method of solving the eigenvalue 
problem automatically yields the eigenfunction \E'o(0), the partial derivatives of eigenvalues are readily obtained by 
taking the quantum average of the corresponding partial derivatives of the potential Vr{Q.,k^). Let us quote the 
derivative with respect to a ( see equation ( |40| ) of section IV. B. for the derivative with respect to /c^) 

(z(L))_ deo{a,k-) ^^^^^^^^^^^^^ ^^^^ 



L 8a 

Once an eigenvalue is known for a value of k^ , one proceeds to the neighboring value fc^ + Ak^ , by using cq + Ak^ 

as a trial eigenvalue and one proceeds to cover step by step the desired range of fc^ = —k^. The fact that the state 

found with this procedure is really the ground state can be checked from the fact that the wave function has no node. 



2. Method b): The transfer matrix iteration 



The search for the smallest eigenvalue to{a, k^) of the transfer operator T{a, k^) defined by eq. (29), ( pOj ) and (|31 
is done by iteration of the mapping |z„) — T |z„_i). Indeed, \zn) = \zq) = \tn){tn\zQ) is reduced to the 
lowest eigenvalue contribution when N 00. Let us write the above mapping as a linear functional transform: 

Znik,0„)^ sin0n-id0n^lT{0n,0n-l) Zn-i{k,0n-i) (51) 



In order to perform the integral over 0n-i we use the following discretization procedure. We divide the variation 
interval < 0n~i < tt into Ug segments '•^""'"^^ < 0n~i < — . The integral over each segment is done with the 
standard Gauss method involving Ug abscissas and Ug attached weights. The integral over the full 0n-i interval is 
then approximated by a discrete weighted sum over d — Ug na points: 



I /(^n-i) sin0„_id0„_i = sin0, 
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With the above integration procedure each iteration step can be reduced to a hnear mapping in a d dimension Euchdean 
space £d- To each function Zn{9) is attached a vector Z„ with d components (Z„)i — ^/wi sin Oi Zn{9i). The additional 
factor has been introduced in such a way that the Hermitian norm (z„|z„) (defined from our discretized integration 
procedure) coincides with the Euchdian norm Z„ • Z„. The Hnear mapping in £d is written as: Z„ — TZ„_i, where 
the elements of symmetric d x d matrix T are given by: 



n-l- 



It is also easily verified that within our finite sum integration procedure we have: (z„|r(a, A:^)|z„_i) = Z„ • TZ 
In order to give a simple criterion of convergence, it is convenient to require that at each iteration step the vector Z„ 
stay on the £d unit sphere. This is achieved by using the following non linear mapping: 

Z„ = ^^"'^ = (52) 
V Z.,i-i • T^Z„_i 

The criterion of convergence is set up as follows : 

Let us define AZ„ = Z„ — Z„_An where An is an integer, typically of the order of one hundred. Z„ is considered 
to be an appropriate fixed point Z/ if ^/ AZn ■ AZ„ < S. In practice we have taken 5 = lO""*; this choice , together 
with An = 200, leads to a very tight convergence test. 

Knowing Zj it is a straightforward matter to get the groundstate energy eo(a, /c^) and its partial derivative with 
respect to a and k^: 

eo{a,k^) = -jlntoia,k^) ^ -^ln(Z/ - TZ/) 

deo{a,k^) ^ A 9r 

9fc2 bto{a,k^) ^'dk^ ^ 

where the Cos matrix elements are given by: {Cos)^ - — Si,j cos 9i. 

It is also possible to get an explicit expression of the eigenfunction ^a{9) from the knowledge of the fixed point 
vector Z f . The eigenfunction obeys the integral equation: 

*o(^)-^/ Ti9,9i)^o{9i)sm9id9, (54) 

To Jo 

Using our finite sum method to perform the integral over 9i in the r.h.s. of the above equation, 5*0 (^) is obtained 
immediately in terms of the components Zf,i of the fixed point vector: 



^0 . 

In our implementation of the method we have taken rtg = 10 . The use of higher values of ng may lead to overfitting. 

In contrast we can vary more freely the number of sectors Ug. The ground state energy £0(0^1^^) and its partial 
derivatives vary by less than one part per million when Us goes from 2 to 8 with Ug = 10. One has to go down to 
Ug = 8 and = 2 to see a variation of about 10"^. We have verified, using the high precision integration subroutines 
provided by the Mathematica software, that the above wave function satisfies, in typical cases, the eigenvalue integral 
equation ( p4| ) to better than 10~^ with the choice rig = 4 and Ug = 10. 

The vector AZ„ = Z„ — Z„_An introduced to set up the convergence criterion can be used to get the energy ei and 
the eigenfunction ^'1(6') of the first excited state with a fairly good accuracy, say better than 10"'^. In particular the 
variation of the energy gap Ae = ci—cq versus 77 will shed some light upon the crossover phenomenon to be discussed 



in section VIII 



VI. MONTE CARLO SIMULATIONS WITH THE DISCRETIZED RLC MODEL 



The Monte Carlo procedure allows to generate configurations of the discretized RLC with a frequency proportional 
to their Boltzmann weight. A full simulation of the discrete model, incorporating the self avoidance effects, the check 



18 



that only unknotted configurations are kept, and the estimation of the Writhe with the non local Fuller formula was 



developed in |g9| in the case of closed DNA chains. Vologodskii and Marko |30| were the first to adapt this formalism 
to the case of a single open supercoiled chain. In order to facilitate the transition from closed chain to open chain, 
they assume that the chain is confined between two impenetrable walls, with the two free ends sitting on different 
walls. From an experimental point of view, one wall is certainly welcome since in the actual experiments one molecule 
end is achored by biochemical links to a glass plate. The other end is biochemically glued to a magnetic bead. In 
the experiments studied in this paper, the bead radius is one seventh of the molecule contour length L, so that an 
unpenetrable wall looks somewhat inadequate to account for the geometrical obstruction of the magnetic bead, even 
if one neglects the very slow processes where the molecule releases the supercoiling by turning around the bead. 
The authors were, of course, aware of the problem and thus they limited their comparison to experiment to relative 
extension larger than 0.3. This wall effect is visible in the limit of zero supercoiling limit where the M.C. results differs 
significantly from the WLC at forces below .IpN ( in that case finite size effects may also play a role since L/A ~ 10 
to 20). If the two conditions < z > / L > 0.3 and F > O.lpN are satisfied then a reasonable agreement was achieved 
with experiments over a rather broad range of forces and supercoiling. The agreement is particularly satisfactory for 
force extension curve at cr = 0.031 and this has allowed the authors to give an estimate of the twist rigidity C, with 
a 20% error bar . We shall return to this last point in section VII. Vologodskii and Marko have also investigated the 
effect of the reduced ioning strength which governs the DNA effective Coulomb radius, which can be characterized by 
the Debye length \d- Passing from 20 mM to 200 mM of NaCl leads to a Debye length reduction by factor VTO- The 
calculated changes in the force extension curves aX a — 0.01 and a — 0.03 are barely visible for stretching forces in 
the range Q.lpN < F < QApN. This includes the {F,a) domain to be considered in the present paper. In the actual 
experiments, the Debye length corresponds to 30 mM of NaCl. Therefore the Vologodskii and Marko findings 
suggest that the self avoidance effects associated with the finite Coulomb radius play a relative minor role in the data 
to be analysed in the next section. 

They present some evidence for the importance of knotting supression in the particular case F = 0.2 pN and 
a — 0.05. A typical simulated knotted configuration is seen to have the ability to absorb the supercoiling more 
efficiently than the corresponding unknotted one: it leaves the chain with a larger extension. Our analysis of the 
F = 0.2 pN hat curve of Fig |l| wiU concern the range \a\ < 0.016 which lies far away from the borderline point 
considered by the authors. 

In the work presented here, we have performed a much less ambitious Monte Carlo simulation. We kept within the 
RLC model without self-avoidance, and the simulation was used essentially as a guide to validate our model, using 
the local formula for xw together with a discrete version of the chain: we checked that it is free of pathology and able 
to reproduce the experimental findings. It suggested to us that the length cutoff b provided by the elementary link of 
the discretized chain could generate the angular cutoff needed to regularize the continuous RLC model, as it is shown 
in section IV. B. Our first computations within the regularized continuous model ( the method a) of subsection V.B.I 
) were checked on few points by Monte-Carlo simulations and the agreement gave us confidence in our approach. 

Subsequently we gave our preference to the transfer matrix iteration method of subsection V.B.2 which uses exactly 
the same theoretical inputs as the Monte Carlo simulation but leads to accurate results with much less computer time. 
The fully discrete model involving three Euler angles per elementary link can be simulated, but it is more efficient 
for our purpose to make use of the fact that the ^ integrals can be done analytically, and thus to work only with the 
two angles 9 and 4> for each link. Indeed, using our computation of section 3 and integrating over the momentum 
k conjugate to the supercoiling angle x after integrating out the ip angles, the partition function of the supercoiled 
DNA molecule can be written as: 

Zix,F)= fv{9,<f>)c^p(^^P^-^(x~ fds<Pil~cos0)] I (55) 




ksT 2L 

This is the path integral which we have discretized and simulated. The discretization procedure is exactly the one 
described in the previous section. We use N elementary links of length b, and the two discretization rules (Ba,E3) for 



the energy function. The path integral measure is substituted by Yln=id4>ndcos9n- The partition function is thus 
expressed as a 2N dimensional integral, and the corresponding probability measure can be sampled by Monte Carlo. 

We use the standard Metropolis algorithm, where a new configuration of the chain is proposed at each step, the 
corresponding change of energy 6E is computed, and the change is accepted with probability min(l, exp(— JE'/fc^T)). 
The point on which one must be careful is the choice of moves. Clearly a choice of local moves (changing 6'„ , 0„ on one 
link at a time) is a very bad one with which a macroscopic change of the molecule takes very many steps, and it leads 
to very long thermalization times. We have checked that it is in effect useless. We need to implement more sizeable 
moves, but moves which do not change the energy too much. One method which we have found rather effective is the 
following. We have first relaxed the constraint involving 9n, which should not change the extensive properties of the 
chain. Then the moves consist in: 
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1. One picks up one link of the chain, number n. 

2. One rotates the section of the chain j G rt + 1, by a global rotation. The rotation, of angle 7 around an 
axis n, is picked up at random, with a uniform distribution for the choice of n on the unit sphere, and a uniform 
distribution of 7 on an interval [—5^ 5] (a value of 5 of order .5 is generally adequate). 

3. One moves to the next link and iterates the procedure. 

This is one of the types of moves which are used in the standard simulations of supercoilcd DNA |2^] . 

We have simulated mostly chains of 300 links, each link b having a length of one tenth of the persistence length. 
We have checked that the finite size effects can be neglected with respect to the statistical errors for this value. For a 
given value of the supercoiling angle x (or rather of its intensive version 77), we perform a simulation, with a number 
of Monte Carlo steps per link of order 10^. The first third of data is used for thcrmalization, the rest is used for 
measuring the distribution of elongation. 

We show in figure || the results for < z > /L, and the statistical error bars. The computations were done 
at C/A = 1.4, and h/A = .10 for values of the force equal to AlQpN , and the amount of supercoiling given by 
T] = 0.0, 0.46, 0.92, 1.39. We see that they are in rather good agreement with the experiment and the transfer matrix 
results. Note that since we do not have to use the two impenetrable walls trick, our computation is also valid in the 
zero supercoiling limit. 

When one increases the degree of supercoiling the thermalization time becomes prohibitive, and an examination 
of the configurations shows that they start building up some small fluctuating plectoneme-like structures. Clearly in 
order to study this regime one must first include self avoidance, and then incorporate moves which are able to shift 

30| ]. We believe that this type of approach can be pursued 
further in order to get a precise comparison with experiments in the strongly supercoiled regime. It is complementary 
to the more restrictive, but more analytical, study at small supercoiling which we develop here. 

VII. ANALYSIS OF THE EXPERIMENTAL DATA ON DNA EXTENSION VERSUS SUPERCOILING 

IN THE SMALL FORCE REGIME. 



In this section we are going to analyse a limited set of data obtained by the LPS-ENS group on the same single 
DNA molecule. They consist of three extension versus supercoiling curves with values of the stretching forces F = 
0A16pN,0A97 pN and 0.328piV. For these data a point-by-point evaluation of the systematic uncertainties is still 
lacking, so only the statistical uncertainties have been considered in the determination of the ratio ^. A simplified 
version of this analysis has been already given by us in a previous publication . 

We have also performed two cuts upon the data in order to exclude from the analysis regions where the validity of 
RLC model is questionable. The first cut allows to neglect the effect of the plane onto which the DNA is attached. 
The second one allows to neglect the self avoidance of the chain. 

First we exclude values of the relative elongation such that i^il^ < 0.1. Experimentally, one end of the molecule 
is attached to a plane which thus implements a constraint z{s) > 0, for the whole chain. We shall see later on that 
when the reduced supercoiling parameter 77 increases from to few units the probability distribution of 9 develops a 
peak near 9 = n. Specially when < i^^i^ < 0.1 , it means that the RLC model is likely to generate configurations 

with < 0. We think that, for the experimental lengths L under study, the regime where is^^ > 0.1 is such that 
typically, the whole chain is in the half plane z > 0. 

The second cut excludes the too high values of the reduced supercoiling angle by imposing the condition ry < 1.5 
for F — 0A97pN and F = Q.328pN. In a following section we shall present evidence that when F = 0.328pN the 
RLC model generates plectoneme-like configurations above a critical value r^c ~ 1 • Our model ignores self avoiding 
effects which are present under the actual experimental working conditions since the DNA Coulombic potential is 
only partially screened. The experimental plectonemes must have a radius |^ larger than the DNA Coulombic radius. 
In contrast, in the RLC model plectonemes with an arbitrary small radius can be generated. For a fixed variation of 
the supercoiling angle, the creation of a plectoneme, having a given length, absorbs, on the average a smaller fraction 
of the DNA chain, in comparison with a situation where self avoiding effects are involved. Indeed we have found that 



'^The plectoneme' radius is the common radius of the two interwound superhelices 



the plectonemes positions efficiently, as was done in |29 
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when the above constraints are not fulfilled the experimental points have a tendency to fall inside the theoretical hat 
curves, i.e for a given 77, experiment gives a smaller ii^ill. 
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FIG. 2. Empirical determination of the cut off length b and the ratio C/A from the hat curves analysis. In the case 
F — .llGpN, we have plotted, versus the cutoff b, the mean ratio {C/A) obtained by averaging the C/A values predicted by 
the RLC model from each empirical hat curve point. The error bar represents for each b the variance ar which leads to a 
measure of the RLC model ability to reproduce the hat curve data. A remarkable agreement is achieved with b — AAA while 
it becomes very poor for b — .06 A. This is consistent with the RLC model singularity near 6 = 0. 



In order to extract the values of C/A and of the cutoff b from the experimental data, we have used the following 
technique. Using equation ( p^ ) the ratio C/A can be written as a function of the reduced supercoiling angle rj and 
the torque k : 



A 
C 



(56) 



With the help of interpolation techniques, one can invert equation (Hq) in order to get k as a function of 



this way, each " hat " curve point of coordinates (ry. 



L ■ 



ML)) 



In 



is associated with an empirical value of the rigidities ratio 



^, once a choice of the cutoff length b has been made. If the RLC model is to give a good representation of the data 
there must exist a value of b such that the empirical values C/A obtained from all points of all hat curves cluster 
around the correct result for C/A. 
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FIG. 3. The elongation versus reduced supercoiling angle, rj = 94.8 cr, for forces F — .116, .197, .328pA'^, from bottom to 
top. The smaller points are the experimental results, the bigger points on the lowest curve are from Monte Carlo simulations 
and the full lines are the predictions of the RLC model obtained as indicated in the text. 

For F =0.116 pN, we have plotted on Figure || the average ratio {C/A) and the variance cr^ = 
y^{{C/A)^) ~ {{C / A) Y versus the cutoff length b. The average and variance are computed over all the 20 ex- 
perimental points of the hat curve. The point with 77 = is an input, used to compute the dimensionless force 
parameter a. Using the scaling properties of the model, we have reduced the data analysis to a two parameters fit: 
h/A and C/A. Since our model reduces to the WLC model in the 77 = limit it is legitimate to use the persistence 
length A obtained within the WLC model from the analysis of the force versus extension curve . As it is apparent 
on Fig. ^ the RLC model with b = 0.14 A leads to a remarkably good agreement with the data. It corresponds to a 
minimum value of 0.03 of the quantity (^^Jj^-^ which measures the ability of our model to reproduce the data. The best 
cutoff length b value is approximately equal to twice the double helix pitch. This is close to the length resolution AZ, 
which we had to invoke in order to justify the assumption of a cylindrically symmetric rigidity tensor. It is interesting 
also to note that the average value (C /A) varies slowly with b. The variance ct^ increases rapidly if one goes to small 
values of the cut-off length ; this is consistent with the fact that the RLC model becomes singular in the limit 6^0. 

Similar results, somewhat less precise, have been obtained for the two other values of F considered in this section. 
They favour the same value of b/A and give values of {C/A) consistent with the ones obtained for F = .WQpN , with 
variances Ur about three times larger. Performing a weighted average upon the whole set of C /A empirical values 
obtained by the three hat curves, taking b — 0.14 A, one gets the following empirical determination of the ratio of the 
two elastic rigidities involved in the RLC model : 

^ = 1.64 ± 0.04 (57) 
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This result is in agreement with the value given in our previous work ^ ~ 1.68. Since we have used data 

which do not incorporate in a quantitative way the systematic uncertainties, the above number should be considered 
as somewhat preliminary, but as it stands, it constitutes a significative improvement upon the previous empirical 
estimates. 

We give in Figure |^ the theoretical " hat " curves using the ratio ^ = (C/A) obtained from the data at each force 
and the cutoff favored value b = 0.14 A. The agreement with the experimental data is very satisfactory and confirms 
the overall consistency of the procedure. It should be stressed that the rather sharp bend connecting a slow quadratic 
decrease to a steep nearly linear falling down, which is observed for 77 ~ 1 on the F — .33 pN hat curve, is rather well 
reproduced by our model. 

Using the values of the persistence length A obtained from the force versus extension curve measured on the same 
single molecule, one can obtain the following empirical value of the twist rigidity: 

C = 84±10nTO (58) 

The statistical error on C is of the order of 2 nm but the systematic errors are expected to be non negligible and the 
overall uncertainty of 10 nm which accounts for those looks reasonable. 

It is here of some interest to quote a recent empirical determination of the twist rigidity obtained from an analysis 
which does not involve the use of the RLC model or the like. A detailed account of the procedure is to be found in 
reference We shall here only quote the result: C = 86± 10 nm, which is in good agreement with the number given 
by equation (|5^). 

Vologodskii and Marko have given previously an estimate of C obtained from one force extension curve at 
a — 0.031, where the agreement with their Monte Carlo simulations is particularly striking : C — 75± 15nm. This 
result agrees with ours. Unfortunately they were not able to get any value of C, leading to an overall good fit for the 
complete set of force extension curves. 

Moroz and Nelson [ pT|j22| have analyzed the "hat" curves, using a perturbation approach to the RLC model. 

Because of the constraint K = — jk'^ > 3, inherent to their method, they explore the force supercoiling angle 

domain given by > .3pN and a < 0.01. This implies that the ranges of relative extension for a given force are 
very narrow : for the top curve of Fig(||) {F = .33pN): 0.72 < (z)/L < 0.75, the higher the force the narrower is 
the interval and closer to unity is its center. It is clear from Fig that the range of relative extension analyzed in 
the present paper is much wider. In fact there is a rather small overlap of the experimental domains involved in the 
two analysis : the intersection of the two data set amounts only to 20 % of our present data set. Moroz and Nelson 
have incorporated hat curves data associated with relatively high force values: 0.6, 0.8, 1.3, S.OpN. For such forces, 
the RLC model is certainly not valid when ct < and even for cr > in the case of the highest force. They J22] have 
derived from their two parameter fit {A, C) the twist rigidity value C — 109 nm (in a preliminary analysis [^l[ based 
upon a more limited set of data they gave C — 120 nm). These authors did not give the uncertainty associated with 
C and the value of A coming out from their two parameter fit. Because of the limited range of relative extension 
values they can fit, we believe that the error bars on their results should be larger than ours. As a comparison, we 
have used C = 109 nm to compute the theoretical hat curves of Fig. ^. As expected, the quality of the fit becomes 
significantly poorer: the overall X2 is multiplied by 5.3 when one goes from C = 86 nm to C = 109 nm. 

VIII. TWIST, WRITHE AND PLECTONEMES. 

In this section we would like to use our solution of the discretized RLC model to study, at a given force (F=.33 
pN), the variations of the torque, twist and writhe thermal averages versus the supercoiling reduced angle rj. As we 
shall see, there exist two very different regimes. Below a rather sharply defined value of the supercoiling angle rjc, 
the DNA chain behaves nearly as an elastic non-flexible rod. Above rjc the torque becomes nearly independent of 
supercoiling while the writhe grows linearly with it. We shall give arguments suggesting that in the high rj regime the 
supercoiling constraint is satisfied by the creation of plcctoncme-like configurations. 



A. Relation between the torque and the average twist 

To begin we are going to prove that the thermal average of the rod twist ( ) is given in terms of the torque F by 
the classical elasticity formula for a non-flexible rod: 

(r.) = ^ = ^^ (59) 



23 



where C = ksT C \& the usual twist rigidity. It is convenient to introduce the joint probabiHty distribution of a DNA 
chain configuration T'{xi, X2) with Tw = Xi xw = Xi\ it can be written as the double Fourier transform: 

7'(xi, X2) oc Z{x\,Xi) = y exp(-ifci xi - 1^2 X2)Z{ki,k2) 

where Z{ki, k2) is given by the functional integrals 

Z{k,M) - / ^) exp U fc2 xw - ^J!SIll±^^\ Zr{k,) 



Zriki) ^ j V{jp) exp (^kiTy, 



^twist 

knT 



k L 

As in section III B, we perforin explicitly the functional integral upon ip and get Zxiki) = exp(— -^). This leads to 
the factorization property: Z{ki, k2) = ZT{ki)Zw{k2), where the writhe term Zw{k2) is given by the path integral: 

Zw{k) ^ j 'D{e,(l)) exp (^ikxw - ^^J^ 

As a physical consequence, the twist T^, and the writhe xw fluctuate independently in the RLC model, when no 
supercoiling constraint is applied upon the free end of the molecule. 

The thermal average ( ), in presence of the supercoiling constraint x — + xiv, is then given by: 

{Tw) = J dxi dx2 Xi Zt{xi) Zw{X2) S{xi + X2 - x) 



1 f dk I dZx \ ~ / , X 

'i—rr- Zw(k) exp{—ikx) 



Z{x, F) J 27r I dk 

Using the explicit form of Zt, one gets the announced formula and by subtraction the writhe thermal average : 



C dx C ksT 

2Lk (9eo(a,-K2) 

.Xw)=X-{Tw) = ^ (60) 



B. The average writhe in the zero stretching force limit: a comparison with Monte Carlo simulations of 

self-avoiding closed supercoiled DNA chains 

The above formulas offer the opportunity to compare, on a specific point, the RLC model predictions with closed 
DNA chain Monte-Carlo simulations |2^, which incorporate self-avoiding effects. Using the same ratio C/A = 1.5 and 
the same value of b/A = .2 as Vologodsku et al. we have computed limi?=o ( Xw)/x =< Wr/ALk > for cr values 
taken in the range: < — cr < 0.04. The comparison with the numbers taken from reference is displayed in Figure 
^. It is apparent that the RLC model computations, and Monte-Carlo simulations are in rather good agreement: 
when \a\ < 0.02 - the range explored in our previous data analysis - the results diverge by less than 8.5% and the 
difference reaches 10% when IctI — > 0.04. They both agree rather well with the measurements performed in references 
p6[ and |27|. If we forget about the possible finite size effects in the Monte-Carlo simulations where A/L ~ 1/12, the 
small deviation may tentatively be attributed to the self-avoiding effects not incorporated in our computations. 
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FIG. 4. Comparison of Monte Carlo 
< Wrjl^Lk >= limj-^o (xw )/x- 



(•) and the RLC model ( continuous line) results for the reduced average 



A torque and writhe versus supercoiling cross-over : a possible sign for thermal excitation of 

plectoneme-like configurations 



In Figure || we have plotted ( for the case F — .33 pN) the torque in ks T unit, k , together with the ratio of 
writhe to twist xw /Tw as functions of the scaled supercoihng variable rj — These theoretical curves have been 
obtained with the same parameters as the "hat" curves of Figure They show a very rapid change of behaviour, a 
quasi-transition, for rjc ~ 1.0, with the following two very different regimes: 



• Below the twist of the DNA chain increases linearly with the supercoiling 77, as in a non flexible rod with an 
effective twist rigidity Cg// — .82 C . The ratio of writhe to twist stays almost constant at the value 0.2. 

• Above rjc the torque depends weakly on the supercoiling 77, the twist becomes nearly constant while the writhe 
increases linearly with the supercoiling. 

The behaviour in the large supercoiling regime is reminiscent of the mechanical instability leading to the formation 
of plectonemes which is easily observed by manipulating macroscopic elastic rods such as telephone cords. One must 
be careful with this analogy because the plectonemic instability corresponds to a zero temperature limit, while in the 
case of DNA, much of the elasticity comes from entropic effects and thermal fluctuations play a crucial role. Yet we 
would like to point towards the existence, in the large supercoiling regime of the RLC, of excitations which share 
some of the properties of plectonemes. 

We define £p as a set of undeformed plectonemes having the axis of their winding helices arbitrarily oriented with 
respect to the force direction ( z— axis ). We then introduce the angular distribution P{9) of the tangent vector t 
about z— axis, averaged along the plectoneme and upon the set £p. Let us sketch the proof of the symmetry relation: 
P{9) — P{tt — 9). Introducing the angular distribution V{ap) of the plectonemes axis and the tangent vector tp/ect(s) 
relative to a vertical plectoneme, P{9) reads as follows : 



P{9)= I V{ap)dapY [ ds 5 (cos 9 - z ■ R {y , ap) ipied 
JQ ^ Jo 
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Ignoring the plectoneme handle contribution, we can rewrite P{6) in terms of the vertical solenoid (/(-dependent tangent 
vector isoi{(j>, 9o) = i?(z, (j))R{y, 9o)z : 

P{6)= V{ap)dap- / dc^{S{cos9~u{ap)-isoi{(l},9o)) + {(l)^TT-(j),9o^TT-eo)) 

Jo inoTT Jq 

where we have used the idendity A • RB = R~^A ■ B and defined the unit vector u(ap) = R{y, —ap)z. Then we 
compute: u{ap) ■ t5o;(0, ^o) = sin ap cos </> sin + cos ckp cos ^J^d by simple inspection, we arrive to the desired 
relation. Because of the (f> averaging, the proof holds true if one takes an arbitrary axis in the z — plane instead 
of the y axis to rotate the plectoneme. In practice the plectonemes are deformed by the thermal Brownian motion. 
More complex structures can also appear, like branched plectonemes. It looks however reasonable to assume that the 
above symmetry property is not affected, on average, by thermal fluctuations. 

In order to measure the degree of symmetry of the distribution of 9 angle along the chain, we introduce the function 
plecto{rj) defined as follows: 

, , , , /;P(g,7?)P(^-g,y?Mcos0) 

plecto{rj) = (61) 

Jo Pi&^v) d{cos9) 

It is clear that plectoij]) reduces to unity for pure plectonemic configurations and goes to zero in the limit of rectilinear 
chains. Within the RLC model, the probability P{9, t]) is given by the thermal average of the molecular axis angular 
distribution when one runs along the molecular chain. Exploiting the quantum mechanics analogy, it is easily proved 
that P{9, rj) is equal to ^'0(6')'^, the square of the ground state wave function introduced in subsection V.B. This has 
been used to compute the function plecto^rj) which is plotted in figure (|^). The shape of P{9,ri) changes in a very 
characteristic way when ones goes from 77 — to « 4. When < < ?7c ~ 1, Pi9, rj) has a rather narrow peak 
at = with a nearly vanishing tail for 9 > ^. As a consequence, the function plecto{r]) is practically null within 
this interval. A secondary peak at 9 — tt begins to develop when rj > rjc and it reaches about the same height as 
the primary peak at 77 « 4. The building of this two bumps structure is accompanied by an almost linear increase 
of the function plecto{'q) which reaches the value 0.9 near rj — A. This behaviour suggests that thermally deformed 
plectoneme-like configurations are responsible for the sharp increase of the writhe-to-twist ratio and the flattening of 
the curve torque versus supercoiling above the critical value 77c « 1. Finally a useful piece of information about the 
quasi-transition near rj > r/c is the study of the variation of the energy gap between the groundstate and the first 
excited state Ae = ei — eo, obtained by the method given in subsection V.B. 2. The corresponding curve appears 
in Figure H under the label "Gap". Ae is a decreasing function of 77, with a rather sudden fall near 77 « 77c. This 
somewhat technical feature has a nice physical interpretation in terms of the correlation length associated with the 
cos 9 fluctuations, at fixed torque , which is given by A/Ae. The jump of this correlation length around 77c, is a further 
sign of a fast change of physical regime. 
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12 3 4 

FIG. 5. Twist, Writhe and Plectonemes in the RLC Model. We display several curves which indicate that a cross-over 
phenomena is taking place near 77^ ~ L After a linear increase, expected for a non-flexible rod, the torque becomes independent 
of the supercoiling above rjc- Simultaneously, the writhe to twist ratio, which was staying constant around .2, starts a rather steep 
linear increase. Near rjc the function plecto{ri), which measures the forward backward symmetry of the t angular distribution 
relative to F, takes off from to the value 0.9.( plecto = 1 for a plectonemic configuration). All these features suggest that the 
cross-over could be attributed to the creation of plectoneme-like configurations. 



The two regimes of supercoiling displayed on Figure ^ are related, within the quantum mechanical formalism of 
sect. V.B.I, to the double well structure of the potential Vr{a, k^, 6) when — —k^. On one hand, the potential has 
a relatively shallow well at = 0, associated with the stretching potential energy. On the other hand, the regularized 
"writhe" potential produces a deep narrow hole at = tt. Let us call ea, and e^, '^i, the eigenvalues and eigenstates 
associated respectively with the semi-classical states localized in the two potential wells. As k is increasing the energy 
levels Ea and Cb are approaching each other. The level crossing is avoided near Kc by a quantum tunneling through the 
finite height wall between the two wells. This accounts for the peculiar variation of Ae near r]c- The mixing of the two 
levels ^'a and ^'b by quantum tunnelling generates the secondary bump in the probability distribution P{0). In the 
vicinity of the near crossing point, the groundstate energy eg depends almost linearly upon the tunneling amplitude 
tab, which has typically a very rapid variation with k^. Since the ratio t^/k is a linear function of ( see equation 
( [43| ) ), -q should exhibit a very sharp increase with k > Kc « 1.4 as is easily seen by turning Figure^ by 90 degrees. 



IX. FLUCTUATIONS OF THE EXTENSION IN A SUPERCOILED DNA MOLECULE 



In this section we would like to give a brief analysis of preliminary measurements of the supercoiled DNA exten- 
sion fluctuations, performed at a force F = 0.33piV by the experimental group in the ENS. We shall compare the 
experimental results with the RLC model predictions usi ng th e parameters b/A and C /A deduced from the extension 



versus supercoiling curves analysis, performed in section VII. The mean square deviations of the fluctuations of the 



molecule extension along the force direction is given by the second derivative of the free energy !F = — ksT log Z{F) 
with respect to F: 
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i5z') = {z')-{zf^{ksTf^^ (62) 

It is convenient to replace the first derivative ^^gp^ by its expression in terms of the average relative extension using 
equation (^6|). The extension mean square fluctuation can then be written as follows: 

We have verified that the value of {5 z"^) obtained in the discretized RLC model in the limit of no supercoiling, 
X = 0, agrees with the one given by the WLC model to better than 3 10^'^. In the situation where the free end of 
the molecular chain is subjected to a supercoiling constraint the fluctuations will be different depending on whether 
the measurement is performed at fixed torque k, or at fixed supercoiling angle x- In the first case the result is readily 
obtained from equation ( ^ ) : 

(^z^)u=-x(^) "''^°j;r''^ (64) 

In the actual experiment the extension fluctuations are measured at fixed x and an extra term has to be added to the 
above expression due to the fact that the torque is now a function of a and ry: 

{Sz )\^^Li^ a ( + 2.-^-^^j (65) 



Using equation (43), the extra term can be transformed by writing 

) de, 



2j|§-(a, —K^) = - — A/C. One gets finally the mean square extension fluctuations at fixed supercoiling angle: 



[5z^)\,^{5z^)U-L[^]a'-l-[^] (66) 



kBT\ ri f dK 

In the above formula the second term is clearly negative so we expect that ((5z^)|,, < ((5z^)|k. The curves giving 
{S z"^) versus the number of supercoils n — w 50 ry are displayed on Figure ^, ( S thick continuous line 

and ( S dashed line. It appears clearly that <C ((5z^)|k, when \ri\ > rjc ^ 1. This implies a strong 



cancellation between the two terms of formula (66), which is then not suitable for an evaluation of )|^. To get 



the thick solid curve of Figure |6|, we have used the fact that the calculation procedure giving, at reduced force a, 
the relative extension versus the supercoiling angle, -^-^(o;,??), is precise enough to allow the evaluation of the partial 
derivative with respect to a by a three points finite difference formula. 

The difference between the two statistical ensembles can be understood qualitatively from the previous section 
considerations and specially by looking at the curves of Figure]^. In a situation where the torque is fixed at k = Kc ~ 1.4 
the supercoiling angle x is practically unconstrained. It can fluctuate rather freely through thermal excitation of 
plectonemes, which have no effect upon the torque. In contrast when 77 is fixed at a value 77 > 77c « 1 the ensemble 
is more constrained. The torque k is still practically fixed at the critical value Kc- As a consequence, the writhe 
~ X - ficL/C has approximately a fixed value and then the only plectonemes which can be created are 
those having a very limited range of writhe. 
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F=0.33 pN 



0.12 




FIG. 6. Mean square fluctuations of DNA supercoiled molecule extension for F — 0.33 pN. The points are the experimental 
data. The thick continuous line is the theory, for the experimental situation where the supercoiling turn number n — « 50 77 
is fixed. The dotted line is the theoretical result for another situation where the torque would be fixed (the number of turns is 
then to be understood as a thermal average value). 

As is apparent on Figure ||, the experimental points agree reasonably well with the RLC model predictions, taking 
into account the quoted uncertainties, which are exclusively of statistical origin. ( Systematic effects are smaller). We 

notice the quasi-critical jump of -^(a, '7), near ric ~ 50, which is also present in the data. 

X. SUMMARY AND CONCLUSION 

Together with new relevant contributions, this paper gives a detailed account of a work of the authors, which 
appeared previously under a very concise letter form [ 12| . It was a rather successful attempt to describe the entropic 
elasticity of supercoiled single DNA molecule in terms of the thermal fluctuations of an elastic rod. 

As in most works based upon a similar model, the rigidity tensor was assumed to be symmetric under rotations 
about the molecular axis and then fully described by two elastic constants, the bending and twist rigidities. Because 
of the DNA helical structure, one may have expected an axially asymmetric rigidity tensor, involving extra elastic 
constants. We have proved in this paper that axial symmetry breaking contributions to the elastic energy are averaged 
out upon a " coarse graining " involving a length resolution AZ about three times the double helix pitch p. Such a 
value corresponds to the actual experimental resolution. The model developed in the present paper is not expected 
to be realistic at length scale below twice the DNA double helix pitch. 

To implement in the RLC partition function the supercoiling empirical constraint x, we have followed a direct 
procedure, instead of adapting to open molecular DNA chains the formalism used in the study of supercoiled plasmids. 
Imposing a well defined regularity condition upon the Euler angles describing locally the molecular chain, we have been 
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able to decompose the empirical supercoiling angle x as a sum of two line integrals and xw , associated respectively 
with the twist and writhe contributions. The local writhe xw is a line integral taken along regular trajectories of 
the chain unit tangent vector t. They are drawn upon the unit sphere 5*2 pierced by a hole having a small radius e, 
localized at the south pole defined by the direction of the stretching force direction. Such a prescription forbids the 
crossing of a spherical coordinate singularity which would invalidate our derivation of the local writhe formula. When 
it is applied in compliance with the above prescription, the local writhe formula allows correct evaluations of the 
writhe of solenoid and plcctoncmc configurations, having an arbitrarily space orientation. More generally the local 
writhe formula, though not explicitly rotation invariant, is shown to lead, in the small e limit, to rotation invariant 
results. 

The partition function of the RLC having a given supercoiling angle x is written as the convolution product of the 
partition function Zt{T^) of a twisted non flexible rod times the partition function Z\y(xw) of a worm like chain 
having a fixed writhe angle. The Fourier transform Zw{k) is shown to be the analytic continuation to imaginary 
time of the Feynman amplitude describing the quantum evolution of unit charge particle moving upon a sphere upon 
the joint action of an electric dc field and a magnetic monopole having an unquantized charge k. The associated 
Hamiltonian HuLcik) is readily obtained from standard Quantum Mechanics rules. 

As in the quantum magnetic monopole problem, we have found a pathology which is associated with the singular 
behaviour of the HRLC'{k) potential term near the south pole. We have proved that an angular cutoff is generated by 
a discretization of the chain involving an elementary link b ( this result was given without proof in our previous work 
[ p2| ) . We have derived analytically from the discretized elastic energy a regularized version of RLC continuous model. 
The singular writhe potential is multiplied by a regulating function going smoothly to zero when sin < ^. We have 

arrived in this way to a well behaved Hamiltonian H'^j^(-;{k), which provided a precise mathematical definition of RLC 
model used in this paper. It leads to a partition function free of any pathology. In an independent work following 
ours ||l^, Moroz and Nelson proposed a high force perturbation method |^l| where no angular cutoff is introduced 
from the start. To finite order, the divergent perturbation series, hopefully asymptotic, does not see the singularity 
at the south pole, though it is always present in the model. In order to make contact with experiment, the authors 
have to impose upon their expansion parameter a sharp " technical " bound. Its effect is to restrict considerably the 
domain of force and supercoiling where a comparison with experiment is possible. This limits the precision of their 
determination of C. 

Within our regularized RLC model, we have developed methods allowing the computation, to an arbitrary precision, 
of the relative extension versus supercoiling curves, the so called " hat " curves. The Fourier transform involved in the 
partition function is evaluated by the saddle point method, in the limit of a contour length L much larger than the 
persistence length A. It leads to a parametric representation of the hat curves in terms of the torque acting upon the 
molecule's free end. The final theoretical ingredient is the ground state energy of the regulated RLC Hamiltonian. It 
is obtained following two methods: a) an explicit solution of the Schrodinger equation associated with H^j^fj(k), h) 
the iteration of the transfer matrix deduced directly from the discretized elastic energy. Though they are not strictly 
equivalent from a mathematical point, they lead to identical results to few %. 

Our data analysis involves three values of the stretching forces: 

F — 0.116, 0.197, 0.328pA^. For each hat curve point and a fixed the RLC model leads to an empirical value of 
the ratio ^ as function of the measured relative extension and supercoiling angle. If the model provides an adequate 
description of the hat curve, the set of empirical values should cluster around the actual value ^. For each cluster we 
have plotted the mean value (^) and the variance versus The best value — 0.14 corresponds to the minimum 
of the ratio crr/(^), which measures the ability of the model to fit the data. The preferred cutoff length &, is found 
to be about two times the double helix pitch and it is close to the length resolution A/ invoked to suppress axially 
asymmetric elastic energy terms. Our determination of ^ = 1.64 ± 0.04 is obtained from a weighted average of (^) 
relative to the three forces involved in the fit, taking ^ = 0.14. This number turns out to be remarkably stable under 
variations of b within the range 0.08A < b < 0.2A: the deviations of ^ from 1.64 stay below the level of 8%. In 
contrast, the quality of the hat curves fit which is very satisfactory at the best b value becomes poorer away from it, 
specially for low b values. Our central value for the twist rigidity C = 84 ± 10 nm differs by about 25% from the ones 
given by Moroz and Nelson The absence of any uncertainty estimate and an analyzed data set having a small 
overlap with ours makes the physical significance of this apparent discrepancy difficult to assess. 

Although these quantities are not yet directly accessible to experiment, we have computed, as functions of the 
supercoiling x, the torque F and the average twist (T^u) (or equivalently the average writhe (xw) — X~ {Tw))- Use 
has been made of a remarkable property: (Tw) is related to F by the linear elasticity formula for a non-flexible twisted 
rod. 

In the case of the highest force {F = 0.33 pN), the two curves, torque and writhe to twist ratio versus supercoiling, 
exhibit a rather sharp change of regime near the same critical value Xc of the supercoiling angle: the torque, after a 
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nearly linear increase, becomes almost supercoiling independent while the writhe to twist ratio, initialy confined to the 
20% level develops a fast linear increase. This behaviour is reminiscent of the buckling instability of a twisted rubber 
tube associated with the creation of plectonemes, able to absorb supercoiling at constant torque. We have shown that 
the configurations excited in the RLC model above the critical value Xc share a simple global symmetry property 
with a set of undeformed plectonemes arbitrarily oriented with respect to the force direction z axis: the angular 
distribution of the tangent vector t has a forward backward symmetry with respect to the z axis. The existence of 
this rather sharp cross over has been confirmed by an analysis of the extension fluctuations versus supercoiling: the 
predicted fluctuations jump near Xc is clearly seen on the preliminary experimental data. 

The overall good agreement of our predictions with the analyzed experimental data seem to indicate that the 
self-avoiding effects, not included in our RLC model, play a limited role in the low supercoiling regime \u\ < 0.02. 
This was suggested by the Monte Carlo simulations of Marko and Vologodskii |^ who use two infinite impenetrable 
walls to adapt the closed chain formalism to the supercoiled open chains. A further test follows from computing 
lim^^o {xw )/x for a values taken in the range: < —a < 0.04. This ratio compares well to the closed DNA chain 
writhe average < Wr/ALk > obtained by Monte Carlo simulations |^^. The two results diverge by less than 10% 
when \a\ < 0.04 and both agree with the experimental data within errors |26|, j2^. It should be said that the two 
calculations use different formulas for the writhe: an open chain local version in the present paper, the non-local loop 
geometrical formula in reference |^ . 

In view of the wealth of experimental information, there are strong motivations to extend the validity of the present 
RLC model to a larger domain of the {F, a) plane. Further work has to be pursued in several directions. Non-local 
constraints in the chain tangent vector t space have to be implemented. There are first the empirical geometrical 
constraints associated with the DNA anchoring glass plate and the finite radius of the tracking bead. The self-avoiding 
effects induced by the Coulomb repulsion within the DNA chain have to be studied for the range of ionic strength 
accessible to experiments. The only practical approach to these problems seem, for the moment, the Monte Carlo 
simulations technique. 

Another interesting perspective is to incorporate in the present model the double helix structure the DNA, in order 
to describe the DNA denaturation transition induced by negative supercoling above F = 0.5 pN. Some first steps 
in this direction have already been taken in |^,^] . Recently a model coupling the hydrogen-bond opening with the 
untwisting of the double helix has been proposed ||l^ . It allows a unified description of DNA denaturation driven 
by thermal fluctuations or induced by the double helix untwisting, in the case a straight line molecular chain. It 
will be of interest to combine this model with the elastic RLC model in order to study the onset of the denaturation 
transition at moderate stretching forces, say below IpN , where bending fluctuations can no longer be neglected. 
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APPENDIX A: ROTATION MATRIX ALGEBRA 



Let us denote by i?(n, 7) the rotation of angle 7 about the unitary vector n. With this notation, the rotation Tl{s), 
which specifies an arbitrary DNA chain configuration in terms of the three Euler angles, is given by : 

n{s) = R (z, R (y, 9{s)) R (z, ^(s)) (Al) 

Another very useful way of writing TZ{s) follows from the rotation group relation: 

R^ R^^R{n,j)Ri = R{R^^h,-/) (A2) 



The proof of (A2) follows from basic properties of a rotation matrix: first, the rotation axis is the rotation matrix 
eigenvector with unit eigenvalue; we verify that it is indeed the case for R^^h: RR^^n — R^^R{h,j)n = i?j~^n; 
second, the rotations R and R(ri,7) have, by construction, the same eigenvalues: 1, exp(i7), exp(— 17) and hence the 
same rotation angle 7. The new form of TZ{s) is then obtained from the simple manipulations: 
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7^(s) = i?(z,(/.)i?(y,0)i?(z» 

= i? (z, (/) + t/.) (z, V') R (y, 0) (z, -0) 

= i?(z,(/. + 7/;)i?(i?(z,-V')y,0) (A3) 

We are going now to discuss the two angular velocity vectors and r2(s) and T(s) defined by relations valid for arbitrary 
vectors X and Y: 

7^(s)7^-l(s)X = n(s) AX (A4) 
7^~^(s)7^(s)Y = T(s) A Y (A5) 



Applying to equation (A5) the simple identity 7?. (a A b) — (TZa) A (TZh) and taking Y — 7?.^^(s)X, one gets im- 
mediately the relation: r2(s) = Tl{s)Y{s). It facilitates the evaluation of the components of r2(s) upon the moving 
trihedron {ei(s)} since we can write: 

= f2(s)-e,(s) = Y(s)-e°(s) (A6) 

To simplify the writing in the explicit computation of T(s), we introduce the notations : Ri = R{z,(f>{s)) , i?2 ~ 
R (y, 0{s)) , i?3 = i? (z, ip{s)). With some elementary matrix algebra we get: 

7^-l(s)7^(s) = {R2Rz)-^R^^RiR2R?. + R^^R2^R2R:i + i?3i?3 

As an intermediary step, we compute : 

i?-ii?-i(n,7(s))7?(n,7(s))i?,X = R-^ (7(s)n A i?,X) = 7i?-in A X 

This result is nothing but the relation ( [A2| ) applied to an infinitesimal rotation. Using the above results we arrive 
finally to the following expression for T(s) 

Y(s) = 0i? (z, -^/') R (y, -61) i + dR (z, -?A) y + ^/-z (A7) 

It is now convenient to decompose T(s) in a longitudinal T||(s) and a transverse part T_l(s) 

T|| = (cos6'0 + V>)z (A8) 
T_L = i?(z,-V') (-0 sin6'x + ^y) (A9) 

We get immediately the quantities appearing in the RLC elastic energy: 

rjg = cos6'0 + V> (AlO) 

= T\ = nl + nl = sin^ e + 9'^ (Aii) 

Our next step is to compute the cylindrical symmetry breaking term Af2(s) = Vl\ — Introducing the angle 
C(s) — arctan {o / (sinflc/))^ , we can write: Tj^ = — _R(z, —ijj — C) x. A physical interpretation of C(s) is obtained 
by computing the s derivative of the tangent unit vector t(s): 

= A t = 7^(s) (Ti A z) 

as 

-^7_L7^(s) i?(z, -1/; - C) (x A z) 
= r!ii?(z,0)i?(y,0)i?(£,-C)y 

We see that — C(s) plays the role of the Euler angle ip{s) vis-a-vis the Seret-Frenet trihedron so that C(s) is clearly 
connected with the writhe. It is now a simple matter to get the component f2i by writing the series of equalities: 

ni = T(s) • el{s) = -n_L (i?(z, ~i! - C) x) • (i?(z, uJas) x) 
= — f7j_x • (i?(z, -0 + C + ^as) x) = — ri_L cos(V-' + C + ^o^) 

In a similar way we obtain fJi = — ilj^ sin('i/; + C + ^os) and we arrive finally at the following expression for Ar2(s): 

A17(s) = nl-VLl = cos 2(?/' -t- C + ^os) (A12) 
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Introducing the length resolution function P(s) = exp(— ^s^/£^) ( Note the change of notation: i stands for 



A I used in the main body of the paper), we proceed with the computation of the average An(s): 
Af2(s) = J dsiP{si - s)AO(si) = J duexp (^-^"^^ ^^(^ + ui) 

Let us first neglect the variation of within the interval (,s — £, ,s + £) and perform a first order expansion in £ 
of the phase ^{s + ut) + <^{s + u£); this is justified since its variation is expected to be of the order of ;j « 0.2 ( We 
have proved explicitly the thermal average inequality: {ip)/r] < 1/A). In this way AQ(s) is transformed into a Gauss 
integral: 

An{s) = An{s)^= J duexp (^-^'^^^ cos (2ue{i} + ( + wq)) 

= An{s)ex.p(^-2f{^ + C + LOof^ (A13) 

Using the estimates |V'|/'^o| ~ ICI/'^o ~ \x\/{L^o) = W\ and the fact that in the present paper our analysis is restricted 
to values of \a\ < 4 10~^, we can write, introducing the pitch p = 2-k/ojo : 

Af2(s) « An(s) exp ( - 2 ( — )' ) (^14) 



Takings = 3.4 nm, t = lOnm we find: « 37; it means that AO(s)/AO(s) is zero for all practical purposes. Using 

instead I = b = 1 nm would not make any difference. Let us say few words about the term (5Ar2(s) involving the 
variation of ^l±^ . A computation similar to the previous one gives the following result: 

dVLi_^ Alii ( IAttI 



5An{s) = — exp ( -^( — sin2(V' + C + wos) 

OS p \ 2 p J 

The rate of variation of 0_l^, which is the inverse of the curvature radius square, is expected to be of the order 1/A 
so that [YZ^ ^gig ( ~ £/A — 0.2. It follows that (5A$l(,s) is still exceedingly small compared to Ari(s). Futhermore it 
is easily shown that the thermal average derivative '^^^g^ ^ vanishes for s ^ A. 

APPENDIX B: THE SYMMETRIC TOP AND THE RLC MODEL 

In this appendix we shall use for convenience a units system where h = c = ksT = 1 and write 9s = t. The elastic 
energy Erlc (eq. (2) and eq. (3) ) is transformed by an analytic continuation towards the imaginary s axis into —i 
times the action integral: 



A{to,ti) = J^' dt ^J2Cjn] + f cos e{t) 

where Ci = C2 = A, C3 = C and f = F/ {ksT). The time derivatives accounts for the relative change of sign between 
Ciflf and the potentiel energy —f cos9{t). The analytically continued partition function of the RLC model is then 
identified with the Feynmann path integral amplitude: 



{9i,(j)i,tpi,ti\9o,(j)o,tpo,to) = / V {0,(p,ip)exp [i / dt£top{t) 



V {0, (p, ip) exp (^i dtC 



The Lagrangian Ctop{t) = \ Qfi? + / cos d{t) describes the motion of a spherical top with inertia moments Ii = Ci, 

under the action of a static electric field Eq. (The molecular electric moment is given by f /Eq.) In order to compute 
explicitly the hamiltonian Titop let us write Ctopit) in terms of the Euler angles and their derivatives: 

= ^ (.^2 sin^ e + e^')+^{7p + ^ COS ef + / cos e (Bl) 
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One gets immediately the conjugate momentums relative to three Euler angles: = A sin^ 94> + cos9p^, = 
C {tp + (j) cos 9) and pg = AO. The symmetric top hamiltonian is then readily obtained: 

C 



H 



top = 2" {'f sin^ ^ + + _ + ^ cos^)2 _ /cos^ 



To get the hamiltonian operator Htop we apply the standard quantization rules: 

d , 9 2-2 ^ 9 . ^ d 

P4> ^ Pel, = , ^ P4> = , P$^P0 = — sme* — 
locp loip sva.9 o9 o9 

To get the partition function sij^O; </'0i ^o, so) we expand the final and initial states upon eigenfuctions 

of the operators and p^: exp{im(j) + ikij)) where k and m are arbitrary real numbers, in contrast with the real 
symmmetric top case, where they are integers. This difference, which follows from a detailed analysis of the physics 
involved ( sec section II B. for details), is responsible for the singular features of the continuous RLC model. The 
partition function with the notations used in the paper to label the initial and final states reads as follows ; 

Z{9{L),(j){L),'4}{L),L\e{Q),Vt) = j dmdkexp{im(f>{L) + iki}{L)) 

{9{L)\exp(^-LIC{k,m)^\9{0)) (B3) 

where the Hamiltonian Kl{k, m) is given by: 

- 1 9 9 (m-cos6lfc)2 fc2 

lC[k,mj = -— — : — - TTTT sin0 — H ^ h -— - ~ fcosU 

^ ' 2Asin6'96i 69 2As\^d 2C ^ 

The experimental supercoiling constraint is implemented by averaging upon 4'{L) and i^{L) the above partition 
function multiplied by the Dirac function 5{<f){L) + V'(i) ^ x)- A straightforward computation gives: 

Z{x) = j dkexp{ikxmL)\exp(^-L}C{k,k)) \9{0)) 
The hamiltonian operator HRLcik) given by equation (16) is recovered by writing: 

HRLc{k) = j(^lC{k,k)-^k'' 

1 d . ^ d ^ fc2 1-COS6I 

:sm9—-acos9+———^ (B4) 



2sin6'a6» 89 2 l + cos6i 

The merit of this direct derivation is to show clearly that k is the unquantized angular momentum of the Euclid- 
ian symmetric top problem associated with the RLC model. It is the breaking of the Quantum Mechanics usual 
quantization rule which is at the origin of the RLC model pathology. 
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